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ABSTRACT 

A fractional-step splitting scheme breaks the full Navier-Stokes equations 
into explicit and implicit portions amenable to the calculus of variations. 
Beginning with the functional forms of the Poisson and Helmholtz equations, we 
substitute finite expansion series for the dependent variables and derive the 
matrix equations for the unknown expansion coefficients. This method employs a 
new splitting scheme which differs from conventional three-step (non-linear, 
pressure, viscous) schemes. The non-linear step appears in the conventional, 
explicit manner, the difference occurs in the pressure step. Instead of solving for 
the pressure gradient using the non-linear velocity, we add the viscous portion of 
the Navier-Stokes equation from the previous time step to the velocity before 
solving for the pressure gradient. By combining this "predicted" pressure 
gradient with the non-linear velocity in an explicit term, and the Crank- 
Nicholson method for the viscous terms, we develop a Helmholtz equation for the 
final velocity. 


LIST OF SYMBOLS 

a = value of non-homogeneous essential boundary condition 

g = value of non-homogeneous natural boundary condition 

A = cross-sectional area 

dA = differential surface area vector 

Bq = 23/12, coefficient used in Adams-Bashforth method 



Bj =—16/12, coefficient used in Adams-Bashforth method 
B 2 = 5/12, coefficient used in Adams-Bashforth method 
B£E C = coefficients 
(Cij)Ibc = coefficients 

D n = Legendre collocation derivative of order n 

Dij = derivative of expansion polynomial j evaluated at node i 

ex, e y , ez = unit vectors in the direction of the co-ordinate axes 

f = body force vector 

E = total number of elements 

M = moment vector 

F = —X 2V 11 * 1 , also represents a force vector 
J[P] = a functional depending on P 
L = periodic length of domain 
L 2 = square-integrable function 
Lk(x) = Legendre polynomial 

it, jt, kt = the maximum number of nodes in the r, s, t directions 
respectively 

I n = a system of interpolating polynomials of order n 

P e = surface integral for side number p and element number e 

S e = total surface integral for element number e 

[J] = Jacobian matrix for the transformation between co-ordinate systems 
n = surface unit normal vector 
s = surface unit tangent vector 
R = position vector 
p = pressure 
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Pi(x) = infinite sequence of orthogonal functions 
p(x), q(x), w(x) = functions in Sturm-Liouville equation 
p = p/p + iV*V, dynamic pressure 
P n = a system of orthogonal polynomials of degree n 
S = an infinite system of orthogonal polynomials 
t = time 

at = time-step size 

t(n) = force per unit area on a surface with unit normal n 
iii = discrete expansion coefficient for the infinite series 
iii = discrete expansion coefficient for the finite series 

V = velocity vector 

V = velocity after the non-linear step 

V = velocity after the explicit viscous step 

V = velocity after the pressure step 

V = corrected velocity after the explicit viscous step 
V ou t = average outflow velocity correction 

U, V, W = x, y, z velocities respectively 

Wi = weight function, integral of the expansion polynomial over domain 
x, y, z = coordinates in global or physical space 

x r = partial derivative of the global co-ordinate x with respect to the local 
co-ordinate r 

r, s, t = coordinates in local or transformed space 

r x = partial derivative of the local co-ordinate r with respect to the global 
co-ordinate x 
Greek and other Symbols 
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a = constant in homogeneous boundary condition 
/? = constant in homogeneous boundary condition 
7 = coefficient in expansion polynomial relationships 
6 = differential of a quantity 
#ij = Kronecker delta 

^i(x) = a combination of Legendre polynomials 

e = infinitesimal quantity 

eijk = alternating tensor 

Q = domain under consideration 

dO = boundary of domain under consideration 

n = OUdO 

<t> = V 11 * 1 

A2= 2/(i*t) 

A = eigenvalue in Sturm-Liouville equation 

A = matrix of coefficients representing A in the expression Ax=b 

II = matrix of coefficients representing b in the expression Ax=b 

p = fluid density 

£ = stress tensor 

a = element surface area 

tTij = i, j component of the stress tensor 

X — moment arm 

V = element volume 

p, = fluid viscosity 

v = fluid kinematic viscosity 

c = Cx e x + Cyey + Cz6z> vortidty vector 
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Subscripts 

o = initial value 
x = free-stream value 
e = element number 
n = direction normal to the surface 
s = direction tangential to the surface 
ijk = co-ordinate system indices 
in = value at inlet 
o U t = value at outlet 
wa ii = value at wall 

XjyjZ = streamwise, vertical, and spanwise values respectively, may also 
refer to partial derivatives with respect to the global co-ordinates x, y , z 
Superscripts 

n = time step number 

e = element number 

i, 2,3,4, 5 , 6 = element side numbers 

INTRODUCTION 

The new splitting scheme, developed by Wessel (1992), contains variations 
on the original three-step splitting method proposed by Korczak and Patera 
(1986). In the previous scheme, the non-linear, pressure, and viscous terms in the 
incompressible Navier-Stokes equations appear in separate fractional steps. By 
introducing intermediate velocities, solutions of these equations yield, 
consecutively, a velocity field based on the non-linear, the pressure and non- 
linear, and the viscous, pressure, and non-linear terms. The final step producing 
the true velocity field. 
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The velocity field resulting from the non-linear step satisfies no boundary 
conditions nor the incompressibility constraint. This velocity field supplies the 
forcing function for the Poisson equation for pressure after applying the 
divergence operator to the pressure step. The intermediate velocity contained in 
the pressure step must satisfy the divergence free constraint, thus it vanishes 
from the Poisson equation for pressure. Instead of solving a second-order Poisson 
equation for pressure, a first-order equation for pressure gradient is solved using 
methods from the calculus of variations-the velocity field follows directly from 
the pressure step. Inviscid boundary conditions on velocity determine the 
pressure boundary conditions; hence, errors of 0(At) occur near solid boundaries. 
Finally, the viscous step employs a Crank-Nicholson scheme yielding a Helmholtz 
equation with a forcing function determined by the velocity after the pressure 
step. Once again a variational form of the governing second-order differential 
equation reduces the order by one. The velocity must satisfy the full, viscous 
boundary conditions; however, it does not satisfy the incompressibility 
constraint. 

The new method varies slightly from the old. Instead of solving for the 
pressure gradient using the non-linear velocity, we include the viscous term from 
the previous time step. Thus, the forcing function appearing in the Poisson 
equation for pressure contains contributions from both non-linear and viscous 
terms. The resulting pressure gradient is not solved for the velocity after the 
pressure step; instead, it, along with the velocity from the non-linear step, and a 
Crank-Nicholson method for the viscous terms, produce a Helmholtz equation for 
the full velocity. The boundary conditions and solution procedure remain 
identical to the original method. The boon comes from including the viscous 
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terms in the pressure gradient prediction, resulting in a quicker solution of the 
pressure step. 
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CHAPTER I 


SPECTRAL APPROXIMATION 

1.1 Spectral Theory 

The expansion of a function u in terms of an infinite sequence of 
orthogonal functions {pi}, u = UiPi, underlies many numerical methods 
of approximation. The most familiar approximation results apply to periodic 
functions expanded in Fourier series. In this case, the i— th coefficient of the 
expansion decays faster than any inverse power of i for smooth functions with 
periodic derivatives. The rapid decay of the coefficients implies that the Fourier 
series truncated after a few terms represents a good approximation to the 
function. This characteristic refers to the "spectral accuracy" of the Fourier 
method. 

Spectral accuracy for smooth but non-periodic functions occurs with the 
proper choice of expansion functions. Not all orthogonal expansion functions 
provide high accuracy; however, the eigenfunctions of a singular Sturm-Liouville 
operator allow spectral accuracy in the expansion of any smooth function, with 
no a priori restriction on the boundary behavior. 

The expansion in terms of an orthogonal system introduces a linear 
transformation between u and the sequence of its expansion coefficients {iii}, 
called the finite transform of u between physical space and spectral space. Since 
the expansion coefficients depend on all the values of u in physical space, they 
rarely get computed exactly; instead, a finite number of approximate expansion 
coefficients result from using the values of u at a finite number of selected 
points— the nodes. This procedure defines a discrete transform between the set of 
values of u at the nodes and the set of approximate, or discrete coefficients. 
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With a proper choice of nodes and expansion functions, the finite series defined 
by the discrete transform represents the interpolation of u at the nodes. 
Maintaining spectral accuracy when replacing the finite transform with the 
discrete transform allows use of the interpolation series instead of the truncated 
series in approximating functions. 

1.2 Sturm-Liouville Problems 

The importance of Sturm-Liouville problems for spectral methods lies in the 
fact that the spectral approximation of the solution of a differential problem 
often occurs as a finite expansion of eigenfunctions of a suitable Sturm-Liouville 
problem. The general form of the Sturm-Liouville problem satisfies 

-SfctPjj) + qa = Awn in (1 £ (-1,1). (1.2.1) 

The real-valued functions, p(x), q(x), and w(x), must behave properly: p(x) 
must be continuously differentiable, strictly positive in (—1,1) and continuous at 
x=±l; q(x) must be continuous, non-negative and bounded in (—1,1); the 
weight function w(x) must be continuous, non-negative and integrable over 
(—1,1). The Sturm-Liouville problems of interest in spectral. methods allow the 
expansion of an infinitely smooth function in terms of their eigenfunctions while 
guaranteeing spectral accuracy. 

1.3 Orthogonal Systems of Polynomials 

Consider the expansion of a function in terms of a system of orthogonal 
polynomials of degree less than or equal to n, denoted by P n - Assume 
(pk}k=o, l- represents a system of algebraic polynomials (with degree of Pk=h) 
mutually orthogonal over the interval (—1,1) with respect to a weight function 
w(x). The orthogonality condition requires 

f Pk(x)p n (x)w(x)dx = 0 whenever m ^ k. (1.3.2) 
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The formal series of a square-integrable function, ueL 2 (— 1,1) J, in terms of the 
system {pk} appears as 

Su = k ? Q iikPk(x), (1.3.3) 

when the expansion coefficients Uk satisfy 

iik = ||pj|| 2 / |u(x)p k (x)w(x)dx. (1.3.4) 

For an integer n>0, the truncated series of u of order n appears as 

P n u = EukPk(x). (1.3.5) 

k = 0 

1.4 Gauss-Lobatto Quadratures and Discrete Polynomial Transforms 

Expanding any u(x)eL 2 (— 1,1) in terms of the coefficients ify, called the 
continuous expansion , depends on the known function u(x). With u(x) is not 
known a priori, a discrete expansion for u(x)— which depends on the values at the 
nodes— must suffice. 

A close relation exists between orthogonal polynomials and Gauss-Lobatto 
integration formulas on the interval [—1,1]. Let Xo,...,x n equal the roots of the 
(n+1)— th orthogonal polynomial p n *i and let w 0 ,...,w n equal the solution of 
the linear system given by 

S o (xj) k Wj = J x^^dx 0 < k < n, (1.4.1) 

where w(x) equals the weight function associated with the Sturm-Liouville 
problem, Wj>0 for j = 0,...,n, and 

.| o p(xj)wj = / jp(x)w(x)dx, (1.4.2) 

hold for all p€P 2 n*i- The positive numbers wj are called "weights" (see Canuto 
et. al.(1988) for proof). This version of Gauss integration produces roots, 


^Identifying the function u(x) as "square integrable" on the given domain 
requires /ju(x)j 2 dx < od. 
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corresponding to the collocation points, which appear in the interior of ( 1,1)- 
Since boundary conditions require one or both end points, a generalized Gauss 


integration formula must include these points. 

The Gauss-Lobatto formula considers 

q(x) = Pn*l( x ) + a Pn(x) + t>Pn-l( x )> (1.4.3) 

with a and b chosen so that q(-l)=q(l)=0. For a given weight function 

w(x) and corresponding sequence of orthogonal polynomials pk, k=0, 1> 2,... 

we denote by x 0 ,...,x n the nodes of the n+1 point integration formula of 
Gauss-Lobatto type, and by wo,...,w n the corresponding weights. 

In a collocation method the fundamental representation of a smooth 
function u on ( — 1,1) appear in terms of its values at the discrete Gauss- 
Lobatto points. Approximate derivatives of the function occur by analytic 
derivatives of the interpolating polynomial. The interpolating polynomial, 
denoted by I n u, belongs in the set P n and satisfies 

I n u(xj) = u(xj) 0 < j < n. (1-4-4) 

Since (1.4.4) represents a polynomial of degree n, it admits an expression given 

by 


I n u = 


n 


S fikPk(x). 
k = 0 


(1.4.5) 


Since the interpolating polynomial must satisfy the function exactly at the nodes, 


we get 

u(xj) = J Q UkPk(xj), (1-4.6) 

where Uk equal the discrete polynomial or expansion coefficients of u. The 
inverse relationship satisfies 

fik = § u(xj)p k (xj)wj, (14.7) 

7k j =0 
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where the coefficients 7k equal 


7 k= f o pi (xj)wj. (1.4.8) 

Equations (1.4.6) and (14.7) enabling transforms between physical space 
{u(xj)} and spectral space {iik} are called the discrete polynomial transforms 
associated with the weights wo,...,w n and the nodes xo,...,x n . 

1.5 Legendre Polynomials 

A collection of the essential features of Legendre polynomials appears 
below. The Legendre polynomials {L k (x), k = 0, 1,...,} equal the eigenfunctions 
of the singular Sturm-Liouville problem given by 

j^l-x^LkM] + *(k + 1)I*M = 0, (1-5.1) 

which equals (1.1.1) with p(x)=l-x 2 , q(x)=0 and w(x)=l. By normalizing 
Lk(x) so that Lk(l)=l, the Legendre polynomials satisfy 

Lk(x) = jkjj gjk (x 2 — l) k (1.5.2) 


and represent the solution to (1.5.1) with boundary conditions (1.5.5). These 
polynomials also satisfy the recurrence relation, expressed as 



Lk*i(x) = xLk(x) fc + 1 Lk-i(x), 

(1.5.3) 

where Lo(x)=l 

and Li(x)=x, 



|L k (x)| < 1, -1 < x < 1, 

(1.5.4) 


L k (*l) = (*l) k , 

(1.5.5) 


|L' k (x)| <|k(k + l), -1 <x<1. 

(1.5.6) 


L'k(*l) = (*l) k |k(k + 1), 

(1.5.7) 

and 

f' L2(x)ix = (k + j)" 1 . 

(1.5.8) 

along with the property that Lk(x) is even if k is even, and odd if k 

is odd. 


The continuous expansion of any ueL 2 (— 1,1) in terms of the Legendre 
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polynomials appears as 


u(x) = E u k L k (x). (1.5.9) 

k = 0 

Multiplying both sides by Lj(x) and integrating from x=- 1 to x=l, gives 

f | u(x)Lj(x)dx = J q / j u k Lj(x)L k (x)dx, (1.5.10) 

where w(x)=l according to (1.3.2). Using the orthogonal properties of the 
Legendre polynomials and (1.5.8) gives 

ij = (j + j) :/ »(x)Lj(x)dx. (1.5.11) 

The Legendre polynomials may appear directly as the expansion functions of 
(1.5.9) or as a combination of Legendre polynomials which satisfy 

^j(xj) = fly, (1.5.12) 


where ip(x) equals 

= n(n+l)L„(xj) x* - S Ln W- (1-5.13) 

This later form allows simpler implementation. The nodes {xj}, j=l,...,n— 1, 
equal the zeros of ^L n (x), with x 0 =- 1, and x n =l. The quadrature weights, 
shown in (1.4.2), satisfy 

J^j(x k )w k = f ^ V'j( x )w(x)dx. (1.5.14) 

Since the Legendre polynomials correspond to w(x)=l, and the expansion 
polynomials satisfy the relation ^j(x k )=5 k j, (1.5.14) reduces to 

X 6 k jw k = f il>\(x)dx (1.5.15) 

k = 0 •'-1 

or 

w j = f i ^j(x)dx. (1.5.16) 

Inserting the expression for the expansion polynomial, (1.5.13), gives after 
integration 


13 



w j “ n(n + l)L"n[5Ej] 2 ’ 0, "’ n (1.5.17) 

for the quadrature weights. The normalization factors 7k, introduced in (1.4.8) 


for the general Sturm-Liouville problem, equal 

7 k = (k -h j)"\ for k < n 
and 7n = 2/n. 


(1.5.18) 

(1.5.19) 


for the specific polynomials given in (1.5.13). 


1.6 Differentiation Using Legendre Polynomials 

Differentiation may occur in either spectral space or physical space. 
Differentiation in spectral space consists of computing the Legendre expansion of 
the derivative of a function in terms of the Legendre expansion of the function 
itself. For example, if u(x)=Ek=o UkLk(x), ^ can be represented as 

^ = = k ? i ^ ik)Lk(x)> ^' 6 ^ 


where 


g^Uk = (2k + 1) p | k+1 iip, P+k odd. (1.6.2) 

The proof of (1.6.2) begins with a relation between Legendre polynomials and 
their derivatives: 

(2k + l)Lk(x) = ^Lk.,(x) - g^Lk-i(x), k>0. (1.6.3) 

Substituting this expression for Lk(x) into (1.6.1) gives 



which upon changing the limits of the summation gives 


(1.6.4) 



- S 

k = -l 


duk* 


+ ^^(X). 


(1.6.5) 


Combining both terms gives 
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as»M = JJi^-iH Z r]E Lk W d (L66) 

where both the terms corresponding to k=0 and —1 vanishes since ^Lo(x)=0 
and jj-L-i(x)= 04 This expression represents the derivative of u(x) in spectral 
space. 

In physical space, the derivative appears as 


m j 

n Q 


as< Su ) = ai u W = k ? 0 4 k ajc Lk(x) 


(1.6.7) 


Equating (1.6.6) and (1.6.7), and recognizing that ^L 0 (x)=0, gives 


k l," k E Lk W _ k ?i 

Since the Lk(x) equal the eigenfunctions of Sturm-Liouville problem, they, along 
with their derivatives, form a linearly independent or orthogonal set. Therefore, 
multiplying (1.6.8) by ^Lj(x) and integrating from -1 to 1 results in 


- duk-i/dx duk*i/dx 

u k = 2k - r 5T+V- 

Evaluating (1.6.2) with n=k+l gives 

(1.6.9) 

j ® 

|^Un-i = (2n - 1) p ? n %, P+n even; 
similarly, (1.6.2) evaluated using n=k— 1 yields 

(1.6.10) 

j ® 

3 ju n *i = (2n + 3) p ? n+2 Up, P+° even. 
Substituting (1.6.10) and (1.6.11) into (1.6.9) gives 

(1.6.11) 

CD ® 

u n = S Up- E Up, p+n even, 

p =n p=n+2 

which reduces to u n = u n . This completes the proof of (1.6.2). 

(1.6.12) 
The two 


derivative expressions given by (1.6.1) and (1.6.7) produce different results in 
practice: 


JFor k=— 1 we have dL-i(x)/dx = c/(l-x J ), where c equals a constant. Since 

the boundary conditions are dLk(±l)/dx=(±l)Mt(k+l)/2 which equals zero for 
k=-l, we see that the constant c must equal zero. 
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^(P n u) i Pn-tg|- 


( 1 . 6 . 13 ) 


The quantity on the left equals the Legendre Galerkin derivative. The error, 
g^(P n u)-P n -ig^ decays spectrally for infinitely smooth solutions. However, for 
functions, u, with finite regularity (not infinitely periodic) this difference decays 
at a slower rate than the truncation error for the derivative — Pn-ig^- Thus 


J UU 

|^P n u) is asymptotically a worse approximation to ^ than P n -q^ (see 
Canuto et. al. (1988)). 

1.7 Legendre Derivatives at the Nodes 

Approximate differentiation in physical space occur by differentiating the 
interpolation I n u (as defined in (1.4.5)) and evaluating it at the nodes. This 
resulting polynomial of degree n— 1, represented as 

Dnti = (1-7-1) 

and called the Legendre collocation derivative of u relative to the chosen set of 
nodes, differs from the Galerkin derivative ^(P n u) since the latter depends on 
the continuous coefficients Uk and the former on the discrete coefficients Uk- 
One method for obtaining the collocation derivative, involves computing 
the values (D n u)(xi), (i = 0,...,n) from the values u(xj), (j = 0,...n), by 
employing (1.4.7) for the discrete Legendre coefficients €Lj, and (1.6.2) for the 
discrete derivative coefficients and computing (D„u)i from 

(D n u)(xj) = k fJj{u k )Vl(xi). (1.7.2) 

A preferred option involves the collocation derivative at the nodes through 
matrix multiplication. It appears as 

(D„u)(x,) = j fix gj(V\(x)] | (1.7.3) 

for i=0,...,n. When DikSgj^(xj), (1.7.3) equals 


du 
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(1.7.4) 


(D n u)(xi) = £ UkDik, for i s= 0,...,n. 
k = 0 

Using (1.5.13) for ^(x) gi ves 

i * k. 

L n (Xk) Xi-Xk’ 

i = k = 0. (1.7.5) 

_(° t Da, i=k = n. 

4 

0, otherwise. 

1.8 Integration Using Legendre Polynomials 

Integration in transform space consists of computing the integral of the 
Legendre expansion of a function. If u(x> E”^ u k ^(x), the integral over the 
domain xe[— 1,1] equals 

J | u(x)dx * f | k E fl iik r{^{x)dx. (1.8.1) 

Assuming the series converges, integration and summation may change places, 
giving 

f u(x)dx « E Uk f ipk{x)dx. (1.8.2) 

Using the integral of the expansion function according to (1.5.16) gives 

J u(x)dx » ^ UkWk- (1.8.3) 
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CHAPTER n 

SPECTRAL ELEMENT METHOD 

2.1 Introduction 

The spectral-element method, a variational procedure in which the 
approximating functions depend on representing the given domain as a collection 
of simple sub-domains, differs from both spectral methods and finite-element 
methods in two ways: (1) pure spectral methods employ high degree 
approximating functions with support defined over the entire domain, and (2) 
finite-element methods use low degree approximating functions with compact 
support (i.e., a given element’s approximating functions differ from zero only 
within the element). Spectral-element methods exploit the advantage of high 
degree functions inherent in pure spectral methods, along with the flexibility 
finite-element methods provide in representing complex domains. The sub- 
domains, or finite elements, equal geometrically simple shapes that permit a 
systematic construction of the approximating functions. These ecumenical 
functions satisfy all boundary conditions and problem data by employing 
concepts of orthogonal polynomials from Sturm-Liouville theory. On an 
demental basis, the dependent variables appear as a finite sequence of the 
approximation functions with coefficients representing the dependent variables at 
a finite number of preselected points (i.e., nodes, whose number and location 
dictates the degree and form of the approximating functions). 

2.2 Partitioning of Domain 

One feature of the spectral-element method distinguishing it from the pure 
spectral method allows representing the given domain by a collection of sub- 
domains. A subsequent transformation maps each sub-domain from the physical 
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(x,y,z) space to the local (r,s,t) space by an isoparametric mapping. The sub- 
domains in local space equal simple geometries, such as cubes in three- 
dimensional space. Two important features in typical geometries dictate this 
mapping: first, the definition of the approximation functions from Sturm- 
Liouville equations only apply to certain well-defined geometries; and second, an 
arbitrary domain cannot accept a collection of simple domains without 
introducing error. By defining the approximating functions element-wise, the 
accuracy of the approximation improves by increasing either the number of 
elements (i.e., refining the mesh) or the degree of the approximating functions. 

In mathematical terms, the total domain d=tt\)dQ splits into a finite 
number, E, of subsets, n e , called finite elements, such that: each fi e is closed 
and non-empty; the boundary dQ e of each f) e is Lipschitz-continuous (no 
singularities, cusps, et cetera); the intersection of any two distinct elements is 
empty, i.e., fi e nflf=®, e^f; and the union H of all elements fie equals the 

total domain, given as 

fi = l fi e . (2-2.1) 

e = 1 

We could not satisfy the last property without the mapping between physical and 
local space. 

2.3 Spectral-Element Interpolation 

By allowing the possibility that each element represents the entire domain 
with the general boundary conditions of the differential equation, the essential 
boundary conditions equal the values of the independent variables at the nodes, 
while the natural boundary conditions get subsumed into the variational form of 
the equation over the element. After assembling the elements, the boundary 
values on portions of the boundaries of elements sharing the boundary of the 
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given domain are replaced by the actual specified values (imposition of boundary 
conditions). 

In the spectral-element method, the minimum degree of the algebraic 
approximating functions depends on the order of the differential equation being 
solved, and the degree of the polynomial in turn dictates the number of 
interpolation points, called nodes, to be identified in the element. 

The approximation functions, also called interpolation functions, depend on 
interpolation of the function and possibly its derivatives at the nodes of the 
element. The nodes placed along the boundary of the element uniquely define the 
element geometry. Place any additional nodes required to define the 
interpolation functions at other points, either in the interior or on the boundary. 
The boundary nodes also enable the connection of adjacent elements by requiring 
equality of the primary degrees of freedom (i.e., variables that appear in essential 
boundary condition) at nodes shared by any two elements. Thus, we cannot 
accurately represent discontinuous primary variables. Such problems arise in, for 
example, the study of compressible flow where shock waves contain velocity 
discontinuities. These functions make poor primary variables in the spectral- 
element model unless we employ special procedures during assembly. 

For each n e , let P£ denote the finite-dimensional spaces spanned by 
linearly independent local interpolation functions {V*?} 1 ? =o of the nodal points. 
Over each element fl e Cfl the approximation I n u e of u e equals 

u e » I n u e = E ut^S(r), (2.3.1) 

1=0 

where the local expansion coefficients u? equal the values of u e at the 
preselected nodes {r?} in the element f} e . As indicated in §1.5, the 
interpolation functions satisfy 
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(2.3.2) 


^i( r j) = 

where <5ij is the Kronecker-delta function. 

2.4 Connectivity (or Assembly) of Elements 

As mentioned earlier, all elements contain boundary nodes defining their 
geometry and allowing connection with their neighbors. The connectivity of 
elements requires equal values of the primary variables in nodes common to 
adjacent elements. Assembling the unique sub-domains into the entire domain, a 
process known as direct stiffness, requires the identification of a universal or 
global system of nodes, and a corresponding set of global expansion coefficients 
for the primary variables. The resulting matrix expression relates the global 
expansion coefficients to the parameters of the governing differential equation 
and boundary conditions. 

As indicated earlier, the primary variables, those associated with the 
essential boundary conditions, appear as the global expansion coefficients of the 
assembled matrix relation. The secondary variables appear in the natural 
boundary conditions. 

2.5 Isoparametric Formulation 

Isoparametric schemes use the same interpolation functions to represent 
both the co-ordinate mapping and the primary variables. Thus, the physical 
space x maps into the local r co-ordinate system by 

x e a I n x e = E xtV^(r), (2.5.1) 

i = 0 

when primary variables appear as 

u e a I n u e = E u?^(r). (2.5.2) 

i = 0 
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CHAPTER HI 

TIME-SPLITTING SCHEME FOR THE NAVIER-STOKES EQUATIONS 


3.1 Governing Equations and Boundary Conditions 

According to the formulation given below, the domain under consideration, 
QcR 3 with boundary dft, may translate but not deform. The Navier-Stokes 
equations on the closed domain, 0=01150, consist of the constant density 
continuity equation given by 

V-V = 0, (3.1.1) 


and the corresponding momentum equation expressed as 

|¥ + (V*V)V = -jE + f+ (3.1.2) 

Introducing 

V*V*V = -(V-V)V + JV(V-V) (3.1.3) 


into the momentum equation gives 

|Y = VxVxV-V(^+iV-V) + f+^V, (3.1.4) 

in which both V and p along with the body force, f, depend on both position 
in the fluid and time. 

The physical boundary conditions for a given problem must allow us to 
divide the entire boundary into regions associated with essential boundary 
conditions (e.g., walls and inlets), natural boundary conditions (e.g., outlets and 
free-streams), and periodicity boundary conditions. The numerical procedure 
handles each of these regions separately. 

The essential boundary conditions appear as 


V = V W all 


(3.1.5) 


on solid wall boundaries, and 


V= V in 


(3.1.6) 
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on inlet boundaries. The natural boundary condition at the outflow equals 

(n-V)V = 0, (3.1.7) 

where n is the outward surface normal vector, and along the free stream 

n-V = n* Vboundary (3.1.8) 

and 

n*(£*n) = 0, (3.1.9) 

where a equals the local stress tensor. Equation (3.1.9) expresses the fact that 
the stress at the free stream must be entirely normal to the boundary. By 
manipulating the pair of conditions (3.1.8) and (3.1.9) we can show that they 
equal the homogeneous natural condition on velocity, (n-V)V=0, (see appendix 
C). The boundary condition along periodic surfaces equals 

V(x+L) = V(x), (3.1.10) 

where L equals the relative position vector between the two periodic boundaries 
All of these conditions refer to velocity; pressure boundary conditions depend on 
the governing equations. Finally, the initial conditions equal V(x,t=0)=Vo(x) 
for xeQcIRa. 

3.2 Splitting Method 

In the variational solution of time-dependent problems, we represent the 
dependent variables in a finite dimensional vector space. The undetermined 
coefficients depend on time, while the base functions depend on spatial co- 
ordinates. This leads to a two-stage approximation, both of which could employ 
variational methods. We choose to discretize the equations in space and iterate 
in time, thus giving a spatial variational problem. Such a procedure, commonly 
known as a semi-discrete approximation, results in a set of ordinary differential 
equations in time. 
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No variational form exists for the full momentum equation; therefore, we 
split it into simple forms and apply variational techniques to each portion 
individually. Employing the splitting scheme followed by Wessel (1992), we 
introduce intermediate velocities, V, V, V, and V, which allows splitting of 
the momentum equation into fractional steps. The scheme employs a "predictor- 
corrector" approach whereby the predicted velocity at time step n+1, which 
results from sequentially computing intermediate velocities based on the non- 
linear terms, the viscous terms, and the pressure terms, determines the pressure 
gradient. The corrected velocity depends on this predicted pressure gradient. A 
brief explanation follows. The first, or non-linear , step appears as 

~ V” = -SPxVxV 11 + P* 1 ; (3.2.1) 

the second, or viscous, step equals 

~ ^°* 1 = i ^ V ”; (3.2.2) 

At 

and the third, or pressure, step includes 

YHH = -V(£ n *V fV^-V 11 * 1 ) = -VP n+1 . (3.2.3) 

The velocity after the pressure step, V 11 * 1 , must satisfy the divergence-free 
constraint. By applying the divergence operator to (3.2.3), a relationship 
between pressure and V" +1 ensues. The solution of this Poisson equation for 
pressure determines the predicted pressure gradient. Using the previously 
computed velocity after the non-linear step and this new pressure gradient, a new 
velocity results after adding the explicit viscous term: 

yn+l ya+l _ _ vpn+ , + ^yjyn (3 .2.4) 

At 

We still must add another $i/V 2 V to the right-hand side to yield the full Navier- 
Stokes equation. We use the Crank-Nicholson method by adding to 
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the velocity V 11 +1 giving an implicit viscous step appearing as 

~ (3.2.5) 

We benefit from this splitting scheme by representing the implicit viscous and 
pressure steps in time-independent, elliptic form: the viscous step in the form of 
H elmh oltz’s equation, and the pressure step as Poisson’s equation. We express 
both equations in variational form and solve them by finding the extremum of 
the corresponding functional. Some of the details of each of the five steps appear 
below. The two steps containing the explicit viscous terms require no exposition. 

3.3 Non-Linear Step 

We solve the non-linear advective term explicitly using a three-step Adams- 

Bashforth method given by (3.3.1) 

V 1141 - = Bo(V*V*V + f) n + B,(VxVxV + f) n ->+ B 2 (VxVxV + f) n ' 2 4 

This hyperbolic operator imposes stability conditions, in the form of a Courant- 
Friedrich-Lewy number, on the scheme. Neither boundary conditions nor 
continuity constraints apply to V 11 * 1 . 

3.4 Pressure Step 

The velocity after the pressure step, V, must satisfy the zero divergence 
constraint. Applying the divergence operator to (3.2.2) gives 

— ) = V-(-VP nt1 )- (3-4.1) 

Since V does not satisfy the zero divergence constraint, (3.4.1) simplifies to 

^‘^ 1>4t = y .VP n+1 . (3.4.2) 

A t 

The boundary conditions on pressure depend on the governing equations. 


JThe three-step Adams— Bashforth coefficients equal Bo=23/12, Bj — — 16/12, 
and B 2=5/12. 
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We obtain the pressure boundary conditions by taking the inner product of 
(3.1.4) with the surface outward unit normal n and rearranging. This yields 
n-VP = -^n-V) + mi-V j V + n-f + n-(VxV*V). (3.4.3) 

This expression represents the exact physical boundary condition on pressure 
obtained from the governing differential equation. Unfortunately, we cannot use 
this form since some of the terms on the right-hand side remain unknown. 
Therefore, the numerical procedure uses a simplified form which neglects the 
viscous term, the body force, and the non-linear term: 

(3.4.4) 

When we discretize this equation in time by writing the time derivative of 
velocity as (V T " l —V n )/bt, we see that unknown terms still remain on the right- 
hand side. By further approximating the time-derivative using intermediate 
velocities we obtain 

n*VP = - V^/at. (3.4.5) 

It appears that we still do not know V 11 * 1 at this point; however, a careful 
vetting of the various boundary conditions reveals otherwise. This mathematical 
approximation to the physical boundary condition may be good or bad depending 
on the characteristics of the flow as indicated below. Where essential boundary 
conditions occur n-V” 41 equals n-V wa n or n*Viniet> and (3.4.5) becomes 

n*VP = — n-(V wa ii — V n+1 )/At. (3.4.6) 

This term differs from the true physical boundary condition at a solid surface (or 
at the inlet if n*(V*V*V)=0) since it lacks both the viscous term, im*V 2 V, and 
the body-force term, n*f. For a flow with no body force, the approximate 
boundary condition differs from the physical boundary condition by the viscous 
term. In other words, the mathematical boundary condition equals the inviscid 
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boundary condition along solid walls. Deville and Orszag (1980) showed that 
this approximation introduces a time-splitting error 0(1) in n-V 2 V over a 
layer of thickness 0{JElv). No error estimate exists at the outlet since even the 
physical boundary conditions remain unknown. 

Along the free-stream the viscous term remains negligibly small under most 
circumstances while the advection term, n-(V*V*V), equals zero, since n 
and V*V*V are perpendicular 4 f Furthermore, when the body force is 
perpendicular to the free-stream boundary— as in most flows— the body-force term 
vanishes, and the simplified boundary condition equals the physical one. 

Along exit boundaries where natural conditions get specified, the terms 
neglected by the approximate pressure boundary condition do not vanish. 
However, their influence remains confined to a region near the boundary , since 
convection sweeps any induced errors out of the domain. 

3.5 Viscous Step 

The implicit viscous step appears as 


JThe viscous term equals n • V 2 V. The velocity can be written in terms of a local 
co-ordinate system where n equals the normal to the surface and s lies 
tangential to the surface. Hence V=VsS+V n n, and n-V 2 V= 
n*[(0 2 V s /ds 2 +d 2 V s /0n 2 )s-l-(5 2 V n /ds 2 +5 2 V n /dn 2 )n]. The first term is zero 
because s and n remain perpendicular. The second term vanishes in most 
common flows. For example, when a homogeneous body force or a body force 
with a vanishing or zero gradient near the free-stream boundary (e.g., the 
"Blasius" body force) forces the flow, the second derivatives of V n vanish 
sufficiently far from disturbance generating structures. Likewise, when the free- 

stream boundary represents a moving plate, as in Couette flow, 0 2 Vn/0S 

usually equals zero, and d 2 V n /0n 2 again vanishes far away from the disturbance 
generating structures. 

t+The boundary conditions specified along a free stream reduce to n-VV=0, 
according to the results of appendix C. Thus, the gradient of velocity lies 
parallel to the surface, or alternatively, the cross product of V and V lies 
parallel to n (the zero normal velocity constraint requires V to lie in the plane 
of the boundary). Therefore, the cross product of V and V*V lies in the plane 
of the surface, or perpendicular to n. Whence n*(V*V*V)=0. 
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yn + 1 _ yn *1 
At 


= W*V** 1 . 


(3.5.1) 


This implicit equation remains unconditionally stable. Therefore, we could avoid 
unreasonable time step restrictions due to the high spatial resolution of spectral 
approximations near the boundaries of elements save the other explicit steps. 

The boundary conditions on the velocity after the viscous step, V“* 1 , equal the 
physical boundary conditions given in §3.1. In formulating the viscous step, we 
did not subsume the zero divergence condition into the expression for V 1141 . 
Hence, V“ 41 does not satisfy the Navier-Stokes equations exactly. Normally, 
this error does not dominate the solution since the velocity after the viscous step 
nearly equals the zero divergence velocity, V”* 1 . Amon (1988) observed that 
the divergence of V 11 41 remains a few orders-of-magnitude smaller than that of 
V 71 * 1 . We transform (3.5.1) into Helmholtz’s equation by adding -V 1141 to 
both sides of (3.5.1) giving 

_Vn*i = _V n * 1 + ji/AtVW 1141 . (3.5.3) 


Rearranging yields 



which appears as 


(3.5.4) 


V2<t»-A2(J> = F, (3.5.5) 

when $=V“ 41 , A 2 =^- and F=-A 2 V n41 . Since $ depends on the velocity 
satisfying the physical boundary conditions, the condition on $ associated with 
essential boundary conditions equals 

(J> = Vwall Of Vjnlet- (3.5.6) 

The exit, or natural, boundary condition equals 

n-V<j> = 0, (3.5.7) 

which also satisfies the free-stream boundary conditions when the free-stream 
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condition on velocity behaves according to the restrictions given in appendix C. 
The periodic boundary conditions equal 4>(x+L)=<f>(x), where L equals the 
relative position vector between the two boundaries. 



CHAPTER IV 


NON-LINEAR STEP 

4.1 Introduction 

The non-linear step, the only explicit part of the three-step time-splitting 
scheme, introduces a time-step size stability restriction. Since only first-order 
derivatives appear, no benefit accrues from casting the governing equation in 
’’weak" form; instead, we apply a collocation variational approach. (When the 
governing equation has an equivalent ’’weak" form, we may multiply it by a 
suitably differentiable test function, integrate over the appropriate domain, and 
then integrate by parts. The resulting variational form contains derivatives of 
lower order than the original equation.) In the collocation technique, the test 
function equals the Dirac-delta function. Since this function has no derivative, 
the resulting solution must contain as many derivatives as the order of the 
governing differential equation. For this first-order equation, therefore, the 
solution occupies the space of functions H^fl). 

4.2 Variational Form 

The Adams-Bashforth three-step procedure applied to the non-linear 
portion of the split Navier-Stokes equation appears in §3.3 and equals 

ynM _yn _ BoAt(V*V*V + f) n + B 1 At(V*V*V + f) n_1 + B 2 At(V*V*V + f) n ‘ 2 . 

(4.2.1) 

To transform this equation into variational form, we multiply by a test function, 
^eH(fl), and integrating over the entire domain, giving 

f [V** 1 -V 11 - BoAt(V*V*V + f) n - BiAt(V*V*V + f) 1 *' 1 - B2At(V*V*V + 

f)”- 2 ]^ = 0. (4.2.2) 

The collocation method uses ipi=6(x.- xi), where the Xi’s equal the locations of 
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the nodal points where the differential equation is satisfied exactly. Introducing 
6(x— Xj) into (4.2.2) and integrating results in the discrete form of the 

differential equation in terms of the unknown expansion coefficients Vyk and 

Vijk- Products, such as V*V*V, also appear in terms of their coefficients at the 
nodes. The result equals 




(4.2.3) 


BoAt(V*V*V + f)? jk + B,nt(V«V«V + f)?ji + B^V-V-V + f)?fi4 
We explore the spatial discritization of the various quantities in the following 
sections. 

4.3 Spatial Discritization of Vorticity 

Express the vorticity, V*V, in Cartesian co-ordinates as 


0 v ,9w dv, , ( d u m^ o ,dv du, 
Vx v= W _ ~fo> ex + _ 


(4.3.1) 


ly ' V 0Z <7X' X (7X tfy 

where the partial derivatives refer to the global co-ordinate system. We must 


transform these derivatives from the global, or physical, co-ordinate system to 


the local (r,s,t) system. According to appendix A, the partial derivatives 
expressed in terms of the local co-ordinates equal 



d e _ T e^ e , .ed e , t edJ 
Jx ~ Tr ft + Sr 3s + tx 3t ’ 

(4.3.2a) 


d e _ T ed e , 5 . t e5 e 

3y + Sy ^s + l m > 

(4.3.2b) 

and 

d e _ T ed e , s ed e t ed e 
Tz + Sz 3s +t£ Bt- 

(4.3.2c) 


(The superscript e refers to the element number.) Introducing these expressions 
into (4.3.1) gives (4.3.3) 


jThe time step number replaces the element superscript for convenience in 
notation; remember, however, that whenever a subscript, such as ijk, represents 
the node number in three dimensions a corresponding superscript should indicate 
the element number under consideration. 
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y v«— f r e^W e , *e^W e , . e ^W e _ e 0V e R e^V e _ t e0VV 

VkV*- (ry^y- + Sy^j- + t y^t Tg FT s *Fs tz FT ^ ei 

, (r edV* + R e^ e 4. t« 5ue r^ 6 - s^ 6 - tggfte 

+( r ^F + Sz cT? +t w ~ Ix JT Sx FT tx FT^y 

, , T e*V« . s e^V e . t e^V* _ r e<?U e _ _ e 0U e _ te 5U e x. 

+( r ^F + Sx ?? +tx Ff r y?F s yjr i m: )ez - 

The three components of the vortidty equal 

r e0W\ 'em*. t e0W e r edV e „0V e _ . e dV e 

Ci= T yjr + s w + t yjr *Z7 Sz W tz Jt ’ 

re dU e , e <5U e , te 5U e r e^W e _ . e 5W e _ t e^W e 
Cy= r W + 8 W + il ~5T - Ix dT Sx Fi~ tr 3T ’ 


and 


- r e5V e , e $V e i t e5V e _ .5U e _ . e 5U e _ e 5U e 

Ci= r *^F + Sx ?s' + tr 3T r ^r s yzr t y?T • 


(4.3.4a) 

(4.3.4b) 

(4.3.4c) 


(4.3.5) 


The x-vorticity evaluated at the local co-ordinate ri,Sj,tk appears as 

(Cx)ijk = (r y )ijk|yijk + ( 8 y)?jk|^ijk + (ty)ijk|^-ijk 

- (r z )?jk|yijk - ( s z)ijk|yijk - (t*)?jk|^ijk- 
The partial derivative with respect to the local co-ordinate, r, equals (see 
appendix A) 

^ijk = Dja^jb^kc- (4.3.6) 

Similar expressions exist for the s and t components. Substituting into (4.3.5) 

g* ves (4.3.7) 

(Cx)ijk = £ [( r y)ijkDia^jb<5kc + (Sy)ijk^iaDjb^kc + (ty)ijk^ia^jbDkc]Wfjk 
abc 

— E [( r z)ijkD ia<5jb^kc + (s z ) ljk^iaDjb^kc + (tz)ijk^ia^jbDkc] Vjjk- 

abc 

Expressing the other vortidty components this way gives an equation for the 
vortidty vector at the local node (ri,Sj,tk): ^4 3 g) 

£ijk = abc{ [l( r y^i k ^ >ia ^ b ^ kc ( s y)>i k ^ iaD j b ^ kc ( t y)1j k ^ia^j bDkc ]Wfj k 

— [( r z)ijkDia< 5 jb< 5 kc + (Sz)ijk^iaDjbftcc + (t z )ijk^ia^jbDkc]Vijk| Cx 
+ |[( r z) ljkDia^jb^kc + (SzJijk^iaDjb^kc + (tz)ijk^ia^jbDkc]Uijk 

— [(rx)ijkDia^jbftcc + (Sj)ijk^iaDjb^kc + (tx)ijk^ia^jbDkc]WijkJ e y 
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+ |j(rx)ijkDia6jtAc + (s x )?jk&aDjb$kc + (tx)ijk^ia<5jbDkc]V1jk 
- [(r y )!jkl>ia^b4c + ( s y)ijk^ia^jb^kc + (ty)tjk^ia^bDkc]UijkjezJ. 

4.4 Spatial Discritization of Cross-Product of Velocity with Vorticity 
The cross-product of the velocity with the vortidty appears as 

V*C= (VCz-WCy)ex + (WCx-UCz)cy + (UCjr-VCx)e*. (4.4.1) 
Using the discrete expressions for the velocity and vorticity gives 

(v-01j k = J c { _ (4 ' 4 ' 2) 

^[(rx)ijkDia^jb^kc + (Sx)ijk<5iaDjb^kc + (tx)1jk^ia^jbDkc]V!jk 

- [(r y )?jkD ia ^jb4c + (Sy)tjk^iaDjb^kc + (ty)?jk«ia«jbDkc)U?jk] V!jk 

- [[(rz)!jkDia^b«kc + (Sz)ijk^iaDjb^kc + (tz)1jk^ia<5jbDkc]Uijk 

- [(rx)ijkDia^jb^c + (Sx)?jk^aDjb^kc + (tx)?jk^ia5jbDkc]W! jk ] W?jk] e x 

+ J^[(r y )?jkDi a 6jb£kc + (Sy)fjk^iaDjb^kc + (ty)ijk^ia^jbDkc]W?jk 

- [(r z )ijkDia<5jb<5kc + (s z )!jk^iaDjb5kc + (tz)ijk^ia^bDkc]V1jkj Wtjk 

- [[(rx)?jkDia«jb^kc + (Sx)?jk^iaDjb4c + (tx)1jk$ia$jbDkc]V?jk 

- [(ry)?jkD ia 5jb4c + (Sy)?jk*iaD jb fcc + (ty)?jk^jbDkc]U?jk]uijk]e y 

•+• ^[(r z )ijkDia^jb<5kc + (s z )ijk<5iaDjb^kc + ( t z) 1 j k ^ ia bD kc] U i j k 

- [(rx)!jkDiaMkc + (Sx)?jk^iaDjb4c + (tx)?jk«ia«jbDkc]W? jk ] U?jk 

- [[(ry)1jkD ia 5jb^kc + (Sy)!jk^aDjb«kc + (ty)?jkMbDkc]W e ajk 

- [(r z )?jkDia^jb4c + (Sz)?jk£iaDjb<§kc + (tz)?jk^ia^jbDkc]V?jkj V!jk] 

for the cross-product of the velocity with the vorticity at the local node (ri,Sj,tk). 

4.5 Summary 

The equation governing the velocity after the pressure step equals 

hli = VJjk + BoAt[(V»C)?jk + %] n ( 4 - 51 ) 

+ B,At[(v-C)?jk + f^jk]"- 1 + B 2 At[(V«C)!jk + *ijk] n ' 2 ; 
where the terms (V*C)ijk appear in (4.4.2). This explicit relation for the 
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unknown velocity coefficients at time step n+1 depends on known quantities 
from the previous step, denoted by the superscript n. Its solution does not 
involve inverting the stiffness matrix; unfortunately, however, this simplicity 
comes with a concomitant loss of accuracy— the collocation scheme searches for 
the solution with a measurement device tolerating first-order errors in time. 
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CHAPTER V 
PRESSURE STEP 


5.1 Variational Approach 

According to §3.4, the governing differential relation for the pressure step 
equals 


y.yn*l 

At 


V-VP. 


(5.1.1) 


This represents an elliptic, Poisson equation; consequently, a corresponding 
variational form exists. For the moment, assume that the functional whose 
Euler-Lagrange equation yields (5.1.1) equals 


J[P] = [_| (VP-VP)At + V n+1 *VP]dV. (5.1.2) 

The standard boundary conditions accompanying a functional of this type equal 

n *|^p = °. (5.1.3) 

where F equals the integrand of the functional. Applying this boundary 
condition formula to the functional shown in (5.1.2) yields 

-n-VPAt + n-V n4l = 0. (5.1.4) 


According to the derivation in §3.4 the appropriate boundary conditions for the 


pressure step equal 

- n-VPAt + n- V 11 * 1 = n-V®* 1 , (5.1.5) 

which does not correspond to that shown in (5.1.4), since an additional n- V n * 
exists on the right-hand side; therefore, we must add a boundary integral to our 
functional. Thus, the functional satisfying the appropriate boundary conditions 
and the governing pressure step relation equals (5.1.6) 

J [P] = Jfj [-} (VP.VP)At + V® 4l *VP]dV - f d(} P(n*V T14l )dcr. 
Applying the Euler-Lagrange equation to this functional yields (5.1.1) with 
boundary conditions (5.1.5). A demonstration of this follows. 
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Up until this point, we accepted the validity of the functional shown in 
(5.1.6) prima facie. We must show that the extremum of this functional does 
indeed yield the differential relation and boundary conditions for the pressure 
step. We begin by considering the variation of J(P] with respect to P: 

*J[P] = Jfj (-VP-V(tfP)At + V“ +1 - V(«>)]dV — (n*V n+1 )fflPdo‘. 

t> • (5.1.7) 

Rearranging gives v 

&J[P] = l-V-(VPSP)it + iPV-VPit + V-(V"* l «P) - «PV-V» ,1 ]dV 

-f m (5-i .8) 

Employing Gauss’ divergence theorem, we transform the first and third volume 
integrals into surface integrals giving 

&J[P] = Jfj [V-VPAt-V-V^aPdV (5.1.9) 

+ f dn [-n-VPnt + n-V^-n-V^^Pda. 

An extremum of J[P] occurs when this variation equals zero. Since the 
variation of P on the boundary of the domain does not depend on its variation 
within the boundary, the extremum occurs when both 

V-VPAt-V-V n * l = 0, (5.1.10) 

within the volume of the domain, and 

— n*VPat + n^V 11 * 1 — n-V 11 * 1 = 0, (5.1.11) 

on the boundary, are satisfied simultaneously. We see that these equations equal 
the governing equation and boundary conditions. Thus, the variation of the 
functional given by (5.1.6) yields the Poisson equation shown in (5.1.1) with 
boundary conditions (5.1.5). 

5.2 Representation in Local Co-ordinates 

We must transform the functional in (5.1.6) from the global x,y,z co- 
ordinate system to the local r,s,t system. The corresponding functional 
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representation for a single element, designated by the superscript e, in local co- 
ordinates equals 

J[p] e = f j f l J J-J (VP • VP)at + V- VP] e | det[J e ] | drdsdt 

- Vj\f\ P c (ii*V e )dAP, (5.2.1) 

where we drop the superscript indicating the time step number and replace it by 
a superscript indicating element number. (Consult appendix B for details of the 
volume element conversion from global to local co-ordinates.) The surface 


integral contains expressions for the summation over the six faces of each 


element. We will consider the details of these terms later. 

This functional, approximated by expressing it in terms of the expansion 

polynomials, for example, P e «E P!bcV’a(r)V'b( s )^c(t), appears as 

abc 

j[p] e « / j c H vtbcPtbc- vibcPibcAt (5.2.2) 

+ VIbc-V|bcPlbc]^ a (r)^(s)^c(t)|det[J e ]|abc<irdsdt -1 P e . 

A brief explanation of the accompanying notation seems appropriate. The 

subscript abc refers to the node r a ,Sb,to while the tilde refers to the 

corresponding unknown expansion coefficients. We postpone treatment of the 

0 

surface integral, conveniently expressed in the interim as IP*, until §5.3. 

Many of the following numerical minutiae receive a more thorough 
treatment in the following chapter dealing with the viscous step. The gradient 
operator in global co-ordinates equals 


VibcPSbc = + 


$Pabc 

~W~ 


ey + 



In the local co-ordinate system this operator equals 


V abc P abc — 


(5.2.3) 


[«||S£(rS), bc + + ^If^tSWlex + 


(5.2.4) 



["'g'f™{ r y)abc + ^g b °( s y)abc + ^^(ty)abc]ey + 

D«be + 4H S ") abc + 

Introducing expansions for the local partial derivatives according to (6.3.5) 
through (6.3.7) gives 

VlbcPIbc = (5.2.5) 

_E^|[(r|)abcDai^bj^ck + (Sx)abc^aiDbj^ck + (tx)abc^ai^bjDck]Cx + 

[( r y)abcDai^bj^ck 4* (Sy)abc^aiDbj^ck + (ty)abc^ai^bj®ck]Cy + 
[( r z)abcDai^bj5ck + (Sz)abc^aiDbj^ck 4- (tz)abc^ai^bj^ck]Czj P?jk* 

The inner product of this term with a term of similar form appears as 

VabcPabc* V|bcPabc = (5.2.6) 

,?k ^[( r x)abcDai^bj5ck + (Sx)abc^aiDbj5ck + (tx)abc^ai^bjl^ck]Cx + 
[( r y)abcDai^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjl^ck]Gy + 
[(rf)abcDai^bj^ck 4* (Sz)abc<$aiDbj^ck + (tz)abc^ai^bjl^ck]Czj Pijk* 
lmn [f( r ^) a ^ c ^ a ^ m ^ cn ^ (Sx)abc^al^bB^cn 4* (tx)abc^al^bnDcn] e x + 
[( r y)abcDal^bm5cn + (Sy)abc^alDbn^cn + (ty)abc^al^bniDcn]Cy + 
[( r z)abcDal^bm^cn + (Sz)abc^all-^biB^cn 4* (tz)abc^al^ba^ cn ]Czl Pfmn- 
Performing the inner product gives 


^[( r x)s 


^abcPabc* VlbcPabc — .? E Pijk Pfnn 
ljk Inn 


(5.2.7) 


rC )abcDai^bj^ck 4* (Sx)abc^aiDbj^ck 4“ (tx)abc^ai^bjl^ck] 
[( r x)abcDal<5bm^cn 4* (Sx)abc^alDbm^cn 4“ (tx)abc^al^bml^cn] + 


l( r y)abcl^ai^bj^ck “4 (Sy)abc^ai^bj^ck 4* (ty)abc^ai^bjl^ck] 
[( r y)abcDal^bm^cn + (Sy)abc^all^biii^cn 4* (tyJabc^al^brnDcu] 4* 
[( r z)abcDai<5bj^ck + (Sz)abc^aiDbj^ck 4* (t£)abc^ai^bjDck] 
[(rz)abcDal^bm^cn "1“ (Sz)abc^alDbm^cn “4 (fz)abc^al^bmDcn]j • 
Introducing these expressions into (5.2.2) gives 
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J[P] e ~ a ^ c V’a( r )V'b( s )^’c(t) I det[J e ] I abc^-i .? k ^Pijk PfmnAt 

^[(rx)abcDai^bj^ck + (Sx)abc^aiDbj ^ck + (tl)abc^ai^bj^ck] 
[(rl)abcDal^bB^cn + (sI)abc^alDbB^cn + (tx)abc$al£b«Dcn ] + 

[(ry)abcDai^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjDck] 
[(ry)abcDal^bm^cn + (Sy)abc^alDbn^cn + (ty)abc^al^b*^cn] + 
[(rl)abcDai^bj^ck + (Sz)ab(AiDbj^ck + (t|) a bc^ai^bjD c k] 
[(r|)abcDal^b»^cn + (Sz)abc<^lDb«^cn + (tt)abc^al<5b»Dcn]j 
+ V|bc* .? k ^[(r|)abcDai^bj^ck + (sl)abc^aiDbj^ck + (t|)abc^ai^bjDck]ex + 

[(ry)abcDai^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjDck]e y + 

[(rl)abcDai^bj^ck + (Sz)abc^aiDbj^ck + (t|)abc^ai^bjD c k]ezj P ijkj dxdsdt 

-IP*. (5.2.8) 

p = l 

Evaluating the integrals by introducing w a = X; ^(r)dr according to (1.5.16) 


gives 


J[P] e » — i £ £ £ P!jk P?«nAtw a wbW c |det[J e ]|abc 

abc i j k Imn 

^[(rl)abcDai^bj^ck + (sl)abc^aiDbj^ck + (tx)abc^ai^bjDck] 
[(rx)abcDal^bm^cn + (Sx)abc^alDbB^cn + (tx)abc^al^bmDcn] + 


[(ry)abcDai^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjDck] 


[(ry)abcDal^bm^cn + (Sy)abc^alDba^cn + (ty)abc^al<5bmDcn] + 


[(rz)abcDai^bj^ck + (Sz)abc^aiDbj^ck + (tf)abc^ai^bjDck] 

[(rl)abcDal^biD^cn + (si)abc^alDb«^cn + (tz)abc^al<5bmD C n]j 

+ £ Vibe* £ P!jkWaWbW c |det[J e ]|abc 
abc ij k 

^[(rl)abcDai^bj^ck + (Sx)abc^aiDbj^ck + (t|)abc^ai^bjD c k]ex + 
[(ry)abcDai^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjDck]e y + 
[(rl)abcDai^bj^ck + (sf)abc^aiDbj^ck + (t|)abc^ai^bjDck]ez| 

_JIPe. (5.2.9) 

p = l 
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Excluding the surface integral, this represents the governing functional in terms 
of the local co-ordinates. 

5.3 Surface Integral 
The surface integral, 

E P* = E f 1 f 1 P e (n-V e )dAP, (5.3.1) 

p = l p-lJ-W-i 

contains contributions from each of the six sides of each element. We must 
express it in terms of the unknown expansion coefficients in the local (r,s,t) co- 
ordinate system. Rearranging (5.3.1) by combining the surface unit normal 
vector, n, with the differential area, dA, gives 

E I pe = E 1 r 1 P e V e -dA p . (5.3.2) 

p = l p=K-K-l 

Introducing these expressions for each of the six faces on the local element gives 

I F e = (5.3.3) 

pxi 

—J jj* j peVMxae* + y3ey + Z3ez)drds 

+ f 1 f 1 P e V e • (x2e x + + Z2ez)drdt 

+f 1 / 1 pe ^ C ‘ ( X3e * + y3Cy + Z3 ^) drds 

-j 1 P e V e -(x 2 e x + y 2 e y + z^)drdt 
—J x j 1 P^- (xie x + yie y + z t ez)dsdt 
p e Ve.( Xl e x + yi e y + ziez)dsdt, 

where we expressed the area vectors shown in (5.3.2) in terms of their 
component parts with the assistance of §A.7 in appendix A. Introducing the 
expansion polynomials into these surface integrals gives 

E P e » (5.3.4) 

p=i 

- f 1 f 1 E P|bi^abr(x 3 e x + y 3 e y + z 3 e z ) a bi^(r)^(s)drds 

+ f 1 f 1 E Pi 2 cVI 2 c-(x2e x + y-fiy + z 2 e z ) a2 cV'a(r)^c(t)drdt 
* -K -l ac 
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+ f 1 f 1 £ Pib^Ibr(x 3 ex + yjey + Z3e*)ab2l&a(r)tf'b(s)drds 

•'-W-l ab 

- f 1 f 1 £ P|icVl lc -(x2e x + y2e y + z 2 e*) a icVfc(r)V'c(t)<irdt 

-K -1 ac 

- f 1 f 1 £ Ptbc^fbc-(xiex + yie y + zie E )ibcV’b(s)lk(t)dsdt 
+ r 1 f 1 S PSbcVSbc-(xiex + yiey + z t e E ) 2 bc^(8)V’c(t)d8dt. 

J -K -1 be 

Evaluating the integrals using w a =/][V’a( r )dr gives 

S P e » (5.3.5) 

p*i 

- £ PlbiVibi-(x 3 ex + y3e y + z 3 ez) a biw a w b 

ab 

+ £ Pa2cVa2c‘ (X2C X + Y’fiy + Z2ez) a 2cW a w c 
ac 

+ £ Pib2VIb2-(x 3 e x + y 3 e y + z 3 e*) a b2w a wb 

ab 

- £ P|lcVIic-(x 2 e x + Tfiy + Z^aicWaWc 
ac 

- £ P1bc^?bc*(xie x + yie y + z^ibcWbWc 

be 

+ £ P$bc^bc-(x,e x + yie y + z 1 e a ) 2 bcWbW c . 
be 

Finally, performing the inner product using VIbc=Uabc€x+V|bc€ y +WIb c ez gives 

EP e » (5.3.6) 

p=i 

- £ PIbi(6lbix 3 + ^lbiy 3 + WIbiz 3 ) ab iw a w b 

ab 

+ £ P|2c(6l2cX2 + ^l2cy 2 + ^S2cZ2)a2cWaW C 
ac 

+ £ P|b2(6lb2X 3 + ^Ib 2 y 3 + W|b2Z 3 )ab2W a W b 
ab 

- £ PIic(fillcX 2 + ^iicy 2 + WllcZ 2 )alcW a Wc 

ac 

- £ Pfbc(6lbcXl + ^Ibcyi + WlbcZl)lbcWbWc 
be 

+ £ P 3 bc(fi 5 bcXi + ^Sbcyi + ^ 3 bcZi) 2 bcWbWc. 
be 
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5.4 Variation of the Functional 

All of the pieces corresponding to the terms in the governing functional 
exist in terms of the local co-ordinates and expansion functions and coefficients. 
Substituting (5.3.6) into (5.2.9) and performing the remaining inner product 
gives 


J[P] e » “i E £ P?jk P?«nAtWaWbW c |det[J e ]|abc 

abcijklan 

^[(rx)abcDaitfbjfick + (Sx)abc^aiDbj^ck + (tx)abc^ai^bj^ck] 
[(rx)abc^al^ba^cn "I" (Sx)abc^alDbm^cn 4 (tx)abc^al^bn^cn] "1“ 


[(^y)abcDai^bj^ck 4* (Sy)abc^ai^bj^ck 4" (ty)abc^ai^bj^ck] 
[( r y)abcDal^bm^cn + (Sy)abc^alDbn^cn 4* (ty)abc^al^bmDcn] 4* 


[( r z)abcDai^bj^ck + (Sz)abc^ai®bj^ck 4* (tl)abc^ai^bj^ck] 

[( r z)abcDal^ba^cn 4" (Sz)abcfcl^bB^cn 4 (tz)abc^al^bmD «]] 

+ E E P?jkw a WbWc|det[J*]|abc (5.4.1) 

abc ijk 

^[(rx)abcDai^bjfick ”1" (Sx)abc^aiDbj^ck 4* (tx)abc^ai^bj^ck]Uabc4 
[(ly)abcDai^bj fick 4" (Sy)abc^ai^bj^ck 4* (ty)abc^ai^bj^ck]Vabc4 


[(r|)abcDai^bj ^ck 4* (Sz)abc^ai^bj^ck 4* (tz)abc^ai^bj^ 



- £ PIbi(U|biX 3 4 Vabiy3 + Wlbiz 3 )abiw a wb 
ab 

+ £ Pa2c(Ua2c x 2 + Vt 2c y2 + Wg 2 cZ2)a2cW a Wc 


+ £ Pab2(Uab2X 3 4 V|b2y3 + W£b 2 Z 3 ) a b2W a Wb 
ab 


- E PIic(U|lcX2 + V|icy 2 + WIi c 2 2 )aicWaWc 

ac 

- E P? bc (UfbcXi + VWl + W!bcZl)lbcWbW c 
be 

+ E P5bc(U5b<oci + VSbcyi + W5bcZi)2bcWbW c . 
be 

The extremum of J[P], or the variation 3J[P]/3P, equals 


42 



dJ[P] e /0P a - E E P!jkAtw a WbW c |det[J e ]|abc 

abc i j k 

^[(ri)abcDai<5bj^ck 4 (sI)abc^aiDbj^ck 4 (tx)abc^ai^bj^ck] 

[( r x)abcDal^bn^cn 4 (Sx)abc^al^bB^cn 4 (tx)abc^al^bB^cn] 4 

[( r y)abcDai$bj^ck 4 (Sy)abc^aiDbj^ck 4* (ty)abc^ai^bjDck] 

[(ry)abcDal^B^cn 4 (Sy)abc^alDbB^cn 4 (ty)abc^al^b»Dcn] 4 

[(ri)abcDai<5bj^ck 4 (s|)abc^ai^bj^ck + (tf)abc^ai^bjDck] 

[(rz)abcDal^bB^cn 4 (s|)abc^alDbB^cn 4 (tz)abc^al^b® Dcn]j 

+ E w a WbW c |det[J e ]|abc (5.4.2) 

abc 

^[(rl)abcDai^bj^ck + (Sx)abc^aiDbj^ck + (tx)abc^ai^bjD c k]UIbc + 

[(ry)abcD a i^bj^ck + (Sy)abc^aiDbj^ck + (ty)abc^i^bjDck]Vtbc + 

[(rl)abcD a i^bj^ck + (Sz)abc^aiDbj ^ck + (t|) a bc^ai^bjDck]W a bcj 

- E (U|bix 3 + V|biy3 + W|biZ3)abiw a wb 
ab 

+ E (Ua2c x 2 4 V|2cy2 4 Wa2c z 2)a2c w a w c 
ac 

+ E (Uab2X3 + V|b2y3 + W|b2Z3)ab2W a Wb 
ab 

— E (U|i c X 2 4 VIicy2 4 Walc z 2)alc w a w c 

ac 

- E (UfbcXi + VWi + W!bcZi)ibcWbW c 

be 

+ E (UBbcXi + V5bcyi + W$bcZl) 2 bcWbW c . 
be 

The middle term in (5.4.2) may appear more conveniently as 

E w a wbw c |det[j e ]|abc (5.4.3) 

abc 

^[(rx)abcDai^bj^ck 4 (Sx)abc^aiDbj^ck 4 (tx)abc^i^bj®ck]U|bc 4 
[( r y)abcDai^bj^ck 4 (Sy)abc^ai^bj^ck 4 (ty)abc^ai^bjDck]Vabc + 
[(rz)abcDai^bj^ck 4 (Sz)abc<£aiDbjfick 4 (ti)abc^ai^bjDck]W|bcj 
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= S Vlbcijk* Vabc^aWbW c | det[J e ] | abc 
abc 

where the divergence operator! equals 

Vabcijk = 

[(rx)abcDai^bj^ck + (Sj)abc^aiDbj ^ck + (tx)abc^ai^bjDck]©x + 
[(ry)abcDai^bj^ck + (Sy)abcfciDbjfick + (ty)abc^ai^bjDck]©y + 
[( r z)abcDai^bj <^ck + (sl)abc5ai^bj^ck + (tS)abc^ai^bjDck]©z- 

Introducing the quantity nf Bna bc> defined by 

Ilfonabc = — | det[J e ] | abc w a w b w c^la^Bb^nC) 

into (5.4.3) gives 


(5.4.4) 


(5.4.5) 

(5.4.6) 


3C) 


£ Vfbcijk ' V|bc w a w b w c | det[J e ] | a bc — — £ £ Vfanijk ’ nfmnabcVfbc 

abc l«n abc 

for the middle term in (5.4.2). The free set of indices in this expression equal 
ijk, while those in the first term in (5.4.2) equal lmn. Rearranging the indices 
gives an equivalent form for the right-hand side of (5.4.6), written as 


£ £ Vijklmn'HijkabcVabc- 

ijk abc 


(5.4.7) 


The extremum of J[P] occurs when the expression in (5.4.2) equals zero. 
Substituting (5.4.7) for the middle term and, setting &J[P]/£P=0 gives 


£ AfmnijkPijkAt — £ £ Vijklan" IlijkabcVabc + S e , 

ij k J ij k abc 

where the surface integral appears as 

S e = 


- E (UIbix 3 + V|biy3 + W|biz 3 )abiw a wb 

ab 


(5.4.8) 

(5.4.9) 


+ £ (Ua2c x 2 + VI 2C y2 + W|2c z 2)a2cW a W c 


ac 


JThe term Vabcijk’ V a bc» obtained from the expression V*V does not equal the 
divergence of velocity. We write the divergence of velocity V- V as VabcijkVijk- 
Hence, the order of indices is important. 
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+ E (U|b2 x 3 + V|b2y3 + W£b2 z 3)ab2 w a w b 

ab 

- E (tSibX 2 + ^Ilcy2 + ^llcZ2)alcWaWc 

ac 

- E (SlbcXl + ^Wl + Wlbc2l)lbcWbWc 
be 

+ E (U5bcXl + VSbcyi + W5bc z l)2bc w b w cj 

be 

and the quantity Atmnijk as 

Af.nijk = WaWbWcl det[J e ] I abc[ (5.4.10) 

[(rl)abcDai^bj^ck + (s£)abcfciD b j^ck + (tx)abc^ai^bjDck] 
[(rx)abcDal^bni^cn + (Sx)abc^alDb«^cn + (tx)abc^al^bniDcn] + 
[(ry)abcDai^bj^ck + (Sy)abc^ail>bj^ck + (ty)abc^ai^bjD c k] 
[(ry)abcDal^bin^cii + (Sy)abc^alDba^cn + (ty) a bc^al^bii)Dcn] + 
[(rl)abc^ai^bj^ck + (s|)abc^ai^bj^ck + (tz)abc^ai^bjBck] 

[( r z)abcDal^bm^cn + (Sz)abc^alDbn^cn + (tz)abc^al^bmDcri]j • 

5.5 Solntion for Pressure 

Equation (5.4.8) applies to a single element. Adding the contributions 

from all elements gives (5.5.1) 

E E z E « 

E E AfnniikPfjkat = E E E V ijklmn* IlijkabcVIbc + ? S . 
e = l ijk J e = l ljk abc e- 1 

In this expression, the pressure coefficients, Pijk. remain unknown; the velocity 


coefficients, V| bc , known from the non-linear step; and the surface integral 

unknown, since it contains unknown velocity coefficients V|bc- By summing 
over all elements, and taking advantage of the continuity of the velocity on the 
faces shared between elements, we greatly simplify the expression for the surface 
integral. Since the surface integrals contain velocities and not derivatives of 
velocity, and, since the velocity remains continuous throughout the domain, the 
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terms cancel on faces shared between two elements. Thus the summation over all 
elements yields a contribution from only those elements on the boundary of the 
domain. Unfortunately, the surface integral contains unknowns which represent 
the velocity after the pressure step. Hence, (5.5.1) contains two unknowns: 

mt W 

pressure Ptjk and velocity Vfjk- Along surfaces where we specify essential 
boundary conditions on velocity, the surface integral depends on known 
quantities; along surfaces where we specify natural boundary conditions, the 

velocity, V?jk, remains unknown. We may circumvent this problem by 
approximating the surface integral of n*V along natural boundaries by n-V, 
since we know V from the non-linear step. However, since V does not satisfy 
boundary conditions, a significant error may accrue from these integrals; 
furthermore, these integrals over the boundaries where both natural and essential 
conditions occur may require extensive computational time. Therefore, we use an 
entirely different technique for the surface integral expression. 

E w 

Since £ S e represents the net flux of velocity leaving the domain, 

e = l 

/ /^q( n* V)dA, we choose to compute the integral of n*V over the inlet and 

exit boundaries and apply the difference of these terms to the exit velocity in the 

form of a correction. Hence, no longer do we satisfy (5.5.1) exactly at each of 

the nodal points. Instead, (5.5.1) holds only in a "global" sense. We 

approximate the surface integral by 

E „ _ 

£ S* Vout* flAout 

( 5 - 62 > 

where V ou t equals the average outflow velocity correction and A ou t equals the 
outflow area. The velocity correction applies to the velocity after the non-linear 
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step according to 


(^ijk)coiTected = V?jk - V ou t, (5.5.3) 

where node ijk and element e correspond to those nodes located along the exit 

portion of the domain, V?jk represents the known velocity after the non-linear 
step, and V 0 ut the average velocity from (5.5.2). We use this corrected 
velocity in (5.5.1) in place of the original uncorrected velocity, and we drop the 

surface integral. The final form equals (5.5.4) 

E _ E z 

X X AfmniikP ijkM = X X X Vijklnn* n?jkabc(Vabc)corrected- 
e = 1 ijk J e=l ijk abc 

Thus, we achieve the standard, linear form, Ax=B, for the unknown pressure 
coefficients. After computing the pressure, the velocity after the pressure step 
results from 

Shjk = Vf jk -V?jkP?jkAt, (5-5.5) 

equivalent to the discritized form of (3.2.2) applied at node ijk of element e. 

Incidentally, a proviso must accompany quantities obtained from 
expressions with derivatives, such as (5.5.5). Since only the variables appearing 
in the essential boundary conditions remain continuous across elements, and not 
their derivatives, the gradient of pressure contains different values at the 
intersection between two elements— one value from each of the two elements. 
Since only one value may exist at any point in the global numbering scheme, 
before equations such as (5.5.5) get solved, we must combine the two values for 
VP at node ijk by taking the average value from the two elements. This seems 

innocuous enough; unfortunately, we break the zero divergence rule on V^jk in 
the process. This may induce important errors in certain problems. 
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CHAPTER VI 


VISCOUS STEP 

6.1 Variational Approach 

The governing differential relation for the viscous step equals Helmholtz’s 
equation, given by 

V^-A^= F, in 0. (6.1.1) 

We employ techniques from calculus of variations to generate a matrix expression 
for the unknown expansion coefficients $ijk- The goal remains finding a 
functional yielding Helmholtz’s equation after application of the Euler-Lagrange 
formula. Assume for the moment, that the functional whose Euler-Lagrange 
equation corresponds to (6.1.1) with homogeneous boundary conditions given by 

o4 + /3(n-V)$ = 0, on dQ, (6.1.3) 

equals 

J[<t>] = f Q [- *grad$:grad$ - iA 2 ^ - $-F]dV. (6.1.2) 

This functional applies to domains with rigid boundaries. In fluid dynamics, 
both homogeneous and non-homogeneous boundary conditions arise. In general, 
the two possible types of non-homogeneous boundary conditions equal Dirichlet 
(essential) when /?=0, given by 

♦ = a, (6.1.4) 

and Neumann (natural) when a=0, given by 

(n.V)* = g. (6.1.5) 

We may incorporate these non-homogeneous boundary conditions in the standard 
functional by adding a boundary integral. 

Since $ must remain fixed along surfaces containing non-homogeneous 
Dirichlet boundary conditions according to (6.1.4), the variation in <|>, written 
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as $*=4+^1 must equal In other words, we impose 

tfdfl) = 0 (6.1.6) 

on surfaces with Dirichlet boundary conditions. We impose no restrictions on rf>, 
along boundaries with natural conditions. Instead, the natural or free boundary, 
where ^ may vary, requires the standard constraint, n*5F/5(grad$)=0, where 
F equals the integrand of the functional. The result, after applying this 
condition to the functional given in (6.i.2), equals n«grad$=0, which does not 
satisfy the non-homogeneous natural boundary conditions. We overcome this 
problem by adding a surface integral to the functional giving 

J[4>] = Jq [ — £grad$:grad$ — f — $*F]dV + 

fffi n4*gdff, (6.1.7) 

where dCl n represents the surface containing natural boundary conditions. The 
variation of the surface integral contributes a term, which when combined with 
n-5F/3(grad$), yields the proper non-homogeneous natural boundary conditions. 
Until now, we accepted the claim that the functional shown in (6.1.7) yields 
Helmholtz’s equation with appropriate boundary conditions without proof. We 
must show that the extremum of this functional does indeed yield Helmholtz’s 
equation and boundary conditions for the viscous step. 

We begin by considering a functional represented as 

Jffl = f n F(x, *, grad«dV. (6.1.8) 

The variation in this functional equals the difference between the functional 
evaluated at 2111(1 4 Thus, <5J=J4*]-J[$] may appear as 

<5J = Jq [F(x, $ + etf, grad$ + cgrad^) 

- F(x, $, grad^)]dV. (6.1.9) 

Expanding the integrand in a Taylor series gives 
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63 = Jn £ *^ + ^gHa$ c:grad ^ dV) 

which, after application of the chain rule, yields 

63 = £ Jn + £ Jn v '(^gFiH$*^ dV - 

Applying this expression to the functional in (6.1.7) results in 

6J = efn [V-grad*- A^-F].*1V + ej V-[-gradf#lV 
+ e fdQ n §****• 

The divergence theorem allows us to express the second integral on the right as 
e f Q V-[-grad$-#lV = - e f dQ n-(grad$-i&)d<7. (6.1.13) 

The integrand on the right-hand side equals zero along the boundary where 
essential conditions apply, since i>=0 in those regions; the only contribution 
comes from boundaries containing natural boundary conditions. Using the 
natural boundary conditions given in (6.1.5) gives 

e Jn V*[-grad$-ifldV = - ej^ n g-^da (6.1.14) 

which exactly cancels the last term in (6.1.12). Thus, the variation in J equals 

63 = c Jn -gratis — A 2 ^ — p ] - (6.1.15) 

At an extremum, we require &J=0 for all admissible ifr, in particular, we 
require &J=0 for all admissible i> which vanish on the surface containing 
essential boundary conditions. Because of the arbitrariness of if inside fi, 

<5J=0 implies 

V 2 $ — A 2 ^ — F = 0, in A. (6.1.16) 

This equals Helmholtz’s equation, confirming our supposition. 

As indicated previously, the variational form incorporates the natural 
boundary conditions in the functional, while the essential boundary conditions do 
not appear. The functional expression makes no reference to a particular domain; 
hence, it applies to both a single element or the entire domain. In the following 


( 6 . 1 . 10 ) 

( 6 . 1 . 11 ) 

( 6 . 1 . 12 ) 
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sections we express this functional in terms of the local co-ordinates on a single 
element. Each element contributes both a volume integral and a corresponding 
surface integral. A functional applicable to the entire domain results from 
s umming over all elements. 

6.2 Representation in Local Co-Ordinates 

The governing functional written in global, or physical space, equals 
J[$] = Jq [- igrad$:grad$ - *A2$-$ - F]dV 

+ Jk i n< N d<r - ( 621) 

We must transform the functional expression from the global (x,y,z) co-ordinate 

system to a local (r,s,t) system, where elements appear as cubes whose local co- 
ordinates range between ±1. The functional representation for a single element in 
local co-ordinates equals 

= J j J j J j[- i grader grad<|> - — 4-F] e | det[J e ] | drdsdt 

+ if\f.r* dAP - (62 ' 2) 

(Consult § A.3 for details of the differential volume conversion between global 
and local co-ordinates.) The surface integral expresses a summation over the six 
element faces. Remember, this term applies only along surfaces characterized by 
natural boundary conditions. Notice that the local co-ordinate system occupies a 
position within an element such that only two co-ordinates vary along a given 
face, while the third co-ordinate remains fixed at *1. The gradient operator in 
the global system equals 

grad$ e = || e x + |y e y + !t e *- (6.2.3) 

In the local co-ordinate system, [r(x,y,z), s(x,y,z), t(x,y,z)], this operator equals 

grad$ e = [|f r £ + ff s * + |f t *l ex (6.2.4) 
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The Jacobian of the transformation equals 

x r x s x t ] e 

J e = y r y s yt (6.2.8) 

Zf Zg 

according to §A.2. Subscripts indicate differentiation with respect to local co- 
ordinates; superscripts define the particular dement under study. 

6.3 Approximating Functions in Spectral Space 

The discrete expansion of a function u(r,s,t) in terms of a finite sequence 
of orthogonal functions of n-th degree approximatdy equals 

u(r,s,t) e a E o ^ J o ufjkA(r)ti( s )fe(t)> (6.3.1) 

where the scalars ufjk equal the discrete expansion coefficients of u(r,s,t) 
rdative to the basis Notice that this expression is not of n-th degree 

unless l=m=n. As indicated in §1.5, differentiation may take place in either 
spectral space or physical space. Normally, we use differentiation in physical 
space. The derivative with respect to the local co-ordinate r in the local space 
equals 

Ur(r,s,t) e a £ E E u!jk^i(r)^j(s)^k(t). (6.3.2) 

1 = 0 J =0 k = 0 
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We do not compute derivatives in this manner since this would entail computing 
analytical expressions for the derivatives of the basis functions; instead, we 
compute derivatives only at the nodal points (r a ,Sb,t c ). This gives 

u r (r a ,3b,t c ) e » 

E # ,| o k £ o U? jk ||i(ra)V’j(Sb)^k(tc). (6-3.3) 

Since the basis functions satisfy both ^i( r a)=^ai and D a i=|^ 1 ( r a)> 
expression reduces to 

Ur(ra>Sb>tc) C 85 ^S^ufjkDai^bj^ckj (6.3.4) 

which may be written more compactly as 

(Ur)Ibc * Uijk Dai^bj^ck- (6.3.5) 

Similarly, derivatives with respect to the s and t directions equal 

(Us)abc # .? k Uijk tfaiDbjtfck (6.3.6) 

and , N 

(ut)abc * ,? k uljk ^ai^bj D ck, (6.3.7) 

respectively. The function u(r,s,t) approximated at the node, (r a ,Sb,tc) equals 

u(r a ,Sb,tc) e « 

E E E u?j k V'i(r a )V'j(sb)V’k( t c) (6-3.8) 

i = 0 j = 0 k*0 

or, in compact form, 

u|bc » E U?jk5ai^bjfick- (6.3.9) 

ljk 

With these expressions for functions and their derivatives, all the ingredients 
necessary to represent the functional in variational form exist. Incidentally, the 
first-order variational form, as opposed to the second-order differential form, does 
not require a second order derivative, clearly an advantage given the complexity 
of the first-order derivatives. 

6.4 Approximation of Products of Functions 
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Since the variational form of our functional contains products of functions, 
some of which are known (e.g., the Jacobian or the body force), while others, 
such as velocity, remain unknown, we need a system to express products of 
functions. The two most common methods of expressing the product of functions 
u and v equal the product of the series, 

u-vs (I n n)(I n v), (6.4.1) 

or the series of the product, 

u-vs In(u-v). (6.4.2) 

The first form yields a polynomial with n J terms since it involves term-by-term 
multiplication of two series; the second method involves a polynomial with n 
terms since the product u and v occurs before the expansion in terms of the 
nodal coefficients. We use the second method. Using the compact notation, the 
product of u and v equals 

u e -v e s .5^ u?jk-v?jkV'i(r)V'j( s )Vv(t). (6.4.3) 

A natural extension of this expression allows a representation of products of three 
or more terms. 

6.5 Surface Integrals 

The term representing the surface integral in (6.2.2) equals 

I e = £11*= £ tfdA*. (6.5.1) 

Using the inhomogeneous boundary condition, g=(n-V)$, gives 

I e = £ fifiA c (n e -V# e dA p , (6.5.2) 

p*i»-i»-i 

which after rearranging gives 

I e = £ f 1 f 1 <k e (n e dA p -V)$ e . (6.5.3) 

p = K - K - 1 

The projection of the differential surface area in the direction of the outward unit 
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normal on each of the six faces equals 

nMA 1 = - 


at t= — 1; 


at s=— 1; 


at t=l; 


at s=l; 


at r=— 1; and 


n e dA 2 = 


n e dA 5 = 


6x 6y Cz 
X r Yr Zr 
L x s y s Z S J 

e* 6y Cg 
X r yr *r 

L x t yt Z t . 


e 

drds, 


drdt, 


Cx 

Xr yr z r 
L x s y s z s J 


e 

drds, 


n e dA 4 = - 


n e dA 5 = 


n e dA 6 = 


Cx Cy eg 
X r yr Z r 

x t yt Z t J 


drdt, 


e x e y eg 
x s y s z 8 
L x t yt z t J 

Cx c y eg 
x s y s z s 
L x t yt z t . 


dsdt, 


dsdt, 


(6.5.4) 


(6.5.5) 


(6.5.6) 


(6.5.7) 


(6.5.8) 


(6.5.9) 


L x t yt ztj 

at r=l. (See §A.7 for details.) Taking the inner product of the projection of 
the differential area in the direction of the outward normal with the gradient 


on face 1; 
on face 2; 
on face 3; 


n e dA‘*V e = — (tx^£ + tyjjy + t^)det[J e ]drds, 

(6.5.10) 

n e dA 2 *V e = — ( s l^ + s JTfy 

(6.5.11) 

n e dA 3 -V e = (tx^ + ty|y + tl^)det[J e ]drds, 

(6.5.12) 

n e dA 4 *V e = ( s *^ + s y^y 

(6.5.13) 
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on face 4; 


n'dA'-V' = - (rSjj + '/ Vf + ^)det[J']ds<lt, (6.5.14) 

on face 5; and 

n e dA 6 *V e = (r^ + r^ + rg^)det[J«]dsdt 1 (6.5.15) 

on face 6, where the derivation of such quantities as t x =dt/dx appear in §A.2. 
Substituting these expressions into (6.5.3) gives for each of the six faces: 


on face 1; 

1“ = - ♦'(«§' + »?$*' + t^*')det[J'ldrds, 

(6.5.16) 

on face 2; 

I!e = + slff)<iet(J']drdt, 

(6.5.17) 

on face 3; 

I" = /.!/.!♦'(«§' + *$' + tsf*)det[J«]drds, 

(6.5.18) 

on face 4; 

!•'= ♦'0*1' + <££ + Gl|| e )det[J e ]drdt t 

(6.5.19) 

on face 5; 

i5e = ♦ e w§' + ^ + rfiVm-kdt, 

and 

(6.5.20) 

on face 6. 

I" = ♦ e (rSf|' + + ^§')det[J']dsdt, 

(6.5.21) 


Substituting the expression for the derivatives 
into the surface integral for face 1, gives 



from (6.2.4) 




+ ** $f)W e det[J e ]drds. 


(6.5.22) 


We must introduce the expansion polynomials and coefficients into these integral 
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expressions. The result for face 1 equals (6.5.23) 



+ t£(r| r+«sr + t| ||)]abidet[J e ]abI+2blV , a(r)V^(s)drds ) 


assuming the quantity within square brackets refers to the nodal points at 
r a ,sb,ti. Similar expressions apply for faces 2 through 6. Substituting 
expressions for the partial derivatives of $ from (6.3.5), (6.3.6), and (6.3.7) 

gives (6.5.24) 

I le S — f 1 f 1 E E [t|(rxD a i5bj£ik + Sx^aiDbj^lk + txtfai^bjDlk) + 

J ~\J • l ab ijk 

ty(ryD a i^bj^ik + S^aiDbj^ik + ty^ai^bjDlk) + tt(rfD a i<5bj^lk + 

sl^aiDbj^lk + t|^ai<5bjDlk)]ab^ijk^Ibldet[J e ]ablV'a( r )^( s )drds. 

The only variables remaining within the integrand equal t^r), V(s), and V<t ); 
therefore, integration reduces to /ii^i(r)dr=wj. Substituting these expressions 

into (6.5.24) gives (6.5.25) 

I le a — E E [tI(rxDai^bj^Xk + s£6 a iDbjtflk + tS^ai^bjDlk) + 

ab ijk 

t y (ryD a i^bj^ik + S^aiDbj^lk + ty^ai^bjDik) + tl(r|D a i<5bj<5lk + 
S^aiDbj^lk + tfMbjD 1 k)]ab4tjk^abldet[J e ] ab iW a W b - 
The superscript tilde identifies unknown discrete expansion coefficients. The 
other terms do not contain the tilde (e.g., tx, r| and det[J e ]) since they 
represent known continuous functions which do not require expression in terms of 
approximate expansion coefficients— these coefficients depend on evaluation of 
functions at the nodes. For completeness, the surface integrals for faces 2 
through 6 equal 

a (6.5.26) 

- E E [s£(rjcD a i£ 2 j£ck + sl^iD 2 jtf C k + tx^i^jDck) + 

ac ijk 

s y (ryD a i£ 2 j£ C k + Sy£ a iD 2 j£ c k + t^aifcjDck) + s|(r|D a i^ 2 j^ck + 
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Sz<5aiD 2 j£ck + t^ai(52jDck)]ab^jk4a2cdet[J e ] a 2cW a w c ; 

I 3e « 

E E [tx(rxDai£bj02k + Sx^aiDbj ^2k + tx5ai^bjD 2 k) + 
ab lj k 

ty(ryD a i^bj 62 k + SyiJai^bj ^2k 4* ty5ai^bjl^2k) 4“ t|(r|D a i^bj^2k + 
Sz^aiDbj^k + tl^ai^jD 2 k)]abl$ijk$ab2det[J e ]ab2WaWb; 

I 4e « 

S # ? [Sx(rSDai5ij5ck + Sx^aiDjj^ck + tx^ai^lj^ck) + 
ac ljk 

Sy(ryD a i^ij^ck *4 Sy^aiDij^ck 4 ty^ai^lj^ck) "4 Sz(rzD a i^ij5ck *4 

Sz^aiD ij(5ck + tz^ai^ljD c k)]abl$ijk$alcdet[J e ]aic w a w c! 

I 5e » 

— E E [r|(r|Da^bj^ck 4 Sx^iiDbj^ck 4* tx^ii^bj^ck) + 

be ij k 

ry(ryDa^bj^ck + Sy^iiDbj^ck + t^aSbjDck) 4" rS(rfDii(5bj5 C k + 

si^iiDbj^ck + t^ii^bjDck)]ib<4ijk$1bcdet[J e ]ibcWbW c ; and 

I 6e ss 

u •?, [^x(^x^2i^bj^ck + Sx^2iDbj4ck + tx&i^bjDck) 

be ij k 

r y (ryD 2 i<5bj <5ck + Sy<5 2 iDbj ^ck + t y ^ 2 i^bjDck) + r®(rfD2i^bj$ck + 

slfeDbjtfck + ti^2i^bjDck)]2b<4iik$!bcdet[J e ]2bcWbWc. 

6.6 Application to Variational Form 

The governing functional approximated in discrete form equals 

J[<M e * f 1 f 1 f 1 S [ — i ^ abc$abc : 7 abc^abc 
w — * labc 

— i^^abc’4*abc — 4*abc* Fabc]V , a( r )V’b(s)V , c(t) | det[J e ] | abcdrdsdt 

+ E P e . 

p=i 

The first term within the integral on the right-hand side appears as 

^abc$abc : ^!bc4abc = 

[(^)abc( r x)ab c + (^)abc(Sx)abc + (^)abc(tx)abc]Cx 


(6.5.27) 

(6.5.28) 

(6.5.29) 

(6.5.30) 


( 6 . 6 . 1 ) 

( 6 . 6 . 2 ) 
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[(|f)abc(ry)!bc + (|f)abc(Sy)!bc + 

[(||)abc(rz)abc + (|f)abc(Sz)abc + (||)abc(tz)abc]ezj * 
[(!f)abc(rx)ibc + (|f)abc(8x)abc + (|f)abc(tx)abc]ex 
[(||)abc(ry)Ibc + (|f)abc(Sy)abc + (|f)abc(ty)Ibc]ey 
[(||)abc(rz)abc + (|f)abc(Sz)abc + (||)abc(tz)abc]eaj • 

Performing the inner product and introducing the discrete expansions for the 
partial derivatives gives 

Vab(4abc:Vab<4abc = .? k jS^ijk'fann (6.6.3) 

^[(rx)IbcDai^bj^ck + (Sx)Ibc^aiDbj^ck + (t x )abc^ai^bjDck] 
[(rx)abcDal^bm^cn + (Sx)abc^alDb«^cn + (t x )abc^l^bii D cii] + 
[(ry)abcDai^bj^ck + (Sy)Ibc^aiDbj^ck + (ty)abc^ai^bjDck] 
[(ry)abcDal^bm^cn + (s y )abc^alDb»^cn + (ty)abc^al^b®Dcn] + 
[(rz)abcDai^bj^ck + (Sz)abc5aiDbj ^ck + (t z )abc^ai^bjDck] 
[(rz)abcDal^bm^cn + (s z )Ibc5alDb«^cn + (t z )abc^al^bmDcn]j • 

Substituting this relation into (6.6.1) gives 

J[<M e * f \f j /j J c V'a(r)t^(s)V , c(t) | det[J e ] I abc[- i.? k jJ n ftk*fcn 

^[(rx)IbcDai^bj^ck + (Sx)abc^aiDbjfick + (tx)abc^ai^bjDck] 

[(rx)IbcDal^ba^cn + (Sx)abc£alDb«£cn + (t x )abc^al^b«Dcn] + 
[(ry)abcDai^bj^ck + (Sy)abc^aiDbj'^ck + (ty)abc^ai^bjDck] 
[(ry)abcDal^bB^cn + (Sy)abc^alDba^cn + (t y )|bc^al^b« D cn] + 
[(rz)abcDai^bj£ck + (s z )abc^aiDbj^ck + (t z )abc^ai^bjD c k] 
[(r z )abcDal^bm^cn + (s z )|bc4lDba^cn + (t z )abc«5al^baiDcn]| 

-iA^bc-4Sbc-^bc-F|bc]drdsdt + p i P®. (6.6.4) 

Again, the only variables within the integrand equal V^r), $ s )> ^CO- 


59 



Using (1.5.16) for /ii^i(x)dx gives 

J[«W e « - i £ .? S ftjk •♦!« | det[J e ] | abcw a w b w c 
a be ij k Inn 

^[( r x)abcDai^bjfick + ( s x)abc^aiDbjfick + (tx)abc^ai^bjDck] 

[(rx)abcDal<5bB^cn 4* (Sx)abc^alDb«ficn 4* (tx)abc^al^bmDcn] 4" 
[( r y)abcDai^bjfick + (Sy)abc^aiDbj^ck + (ty)abc^ai^bjDck] 
[(ry)abcDal^ba^cn + (SyJabc^alDbnficn + (fy)abc^al^b»Dcn] + 
[( r z)abc^ai^bj^ck + (Sz)abc^aiDbj^ck 4- (tz)abc^ai^bj^ck] 
[(rz)abcDal<5bm^cn + (Sz)lbc^all^bm^cn + (tz)abc^al^bm^cn]j 

-WE ft 

mn*ftmn | det[J e ] | lanWlW.Wn 

lmn 


- S ^! Bn -F! Bn |det[J e ]|i Bn wiw.w n + E I pe . (6.6.5) 

lmn p*l 

The equation which minimizes the functional J[$] occurs when 63/ 5$i mn = 0, 
and equals 


£[( r x)abcD 


* “ , 5 c ifk ^ k l det ( J 'll» b ' W »* bWc 


(6.6.6) 


^ai^j^ck + (Sx)abc^aiDbj^ck + (tx)abc^ai^bj^ck] 
[(rx)abcDal^bm^cn + (Sx)abc^alDbm^cn + (tx)abc<5al^bmD cn ] + 


[( r y)abcDai^bj<5ck 4* (®y)abc^ai®bj ^ck + (ty)abc<$ai^bjDck] 
[( r y)abcDal^bm^cn 4* (Sy)abc^al^bB^cn "4 (tyjabc^al^m^cn] + 
[(rz)abcDai^bj^ck 4* (Sz)abc^aiDbj^ck 4“ (tz)abc^ai^bjDck] 
[( r z)abcDal^bm^cn 4* (Sz)abc^alDbm^cn + (tz)abc^al^bmD cn ]j 
-A 1 # mn | detfJ 6 ] | lmnWlW B W n - Pf Bn |det[J e ]| lunWlWnWn + 

JJdn n ‘^glla^ a + lif(p? 1 Ipe ]- 

The last two terms cancel.}: Setting fiJ/% Bn =0 gives a more compact form of 


{The integrand of the second-last term equals n-5F/dgrad$ which also equals 
— n*grad$ for the functional defined in (6.2.2V The integration of this term 
takes place over both the essential and natural boundaries. The variation of the 
surface integral with respect to <j> equals g integrated over the boundary where 
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(6.6.6): 


S [Afomijk - A 2 nf B nijk]#jk =.? n?.nijkF1 ik , (6-6.7) 

ljk 

where II and A appear as 

nfomijk = — |det[J e ] I ijkWiWjWk^i^j^k (6.6.8) 

and A? B nijk = S |det[j e ]|abcw*wbw c f (6.6.9) 

abc l 

[(rx)abcDai^bj^ck + (s x )&bc5ai^bj^ck + (tx)abc^ai^bjDck] 
[(rx)abcDal^bB^cn + (8x)abc^alDb*^cn + (tx)abc^al^b«Dcn] + 
[(ry)abcDai^bj^ck + (Sy)abc^aiDbj fick + (ty)abc^ai^bjDck] 
[(ry)IbcDal^bm^cn + (Sy)Ibc^alDbn^cn + (tyJabc^al^bmDcn] + 
[(r z )abcDai^bj^ck + (Sz)Ibc^aiDbj^ck + (t z )abc^ai<5bjDck] 

[(r z )IbcDal^bm^cn + (s z )abc^alDbB^cn + (tz)atxAl^bnD C n]j • 

Performing the multiplication associated with (6.6.9) gives 

A?mnijk= E det[j e ]|abcw a wbw c [ (6.6.10) 

abc 

(r®r| + Tyry + r|ri)abcDai^bj^ckDal^bB^cn + 

(r*Sx + rySy + rts|)abcDai^bj ^ck^alDba^cn + 

(rltx + rft£ + r|t|)abcDai^bj^ck^al^b«iDcn + 

(s£r| + Syr® + s|r|)abc^aiDbj ^ckDal^b*<5cn + 

(s|Sx + SySy + sisi)abc£aiDbj£ck^alDba£cn + 

(sitx 4" Syty + S^t|)abc^aiDbj <5ck^al^b»Dcn 4" 

(t£r® t tyly + t|r|)abc^ai^bjDckDal^b«^cn 4" 

(txSx "f tySy + t|s|)abc^ai^bjI^ck^alI^bB^cn 4 

(txtx 4" tyty 4" tft|)abc^ai^bi^ck^al^bB^cn ]• 


natural conditions exist. Since g equals n*grad<), these terms cancel along the 
boundary where the natural conditions apply; the variation of ^ along the 
boundaries where essential boundary conditions applv equals zero. Whence, 
these terms contribute nothing to the variation of J[<H- 
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This expression simplifies after contracting on the free indices associated with the 
Kronecker delta; this requires several steps. First, multiply (6.6.10) by $ijk 
and contract on the indices a, b, and c giving 

.? Af.nijkfek= (6.6.11) 

ijk 

£ E(rxrx r»Iy + ^E)sinDsi^aj4ikllal | | un^a^n^n^ijk 

ijk a 

+ £ £(r|s| + rSsy + r|si)ibnDii<5bj^nkDb.|det[J e ]|ibnWiw b w n ^!jk 

ijk b 

4 E E(rxt| 4 Tyty 4 r|tt)lncDli£«j£ckDcn| d€t[J e ] | lBc w l w mW<$ijk 

ijk c 

4 E E(SxrI 4* Syr® 4 s|r|)amn^ai®Bj^nkDal|d6t[J e ] | amn w a w m w n$ijk 

ijk a 

+ £ £(s|s| + sSs$ + S|si)lbn^liDbj^nkDb.|det[J e ]|lbnWlW b Wn^jk 
ijk b 

+ £ £(s|t| + sStS + s|t|)i 1 ncftiD B j5ckD C n|det(J e ]|i.cWiw 1B wc$?jk 
ijk c 

4 E E(txr| + tyly 4 tfr®)aiDn^ai^injI^nkDal|d6t[J e ] | amnWaWmWn^ijk 
ijk a 

+ £ £(tlsi+ tSsS + t|si)ibn.^ii5bjDnkDb«|det[J e ]|ibnWiw b w n ^jk 
ijk b 

+ £ £(t|t| + tyt^ + tftDiBc^li^mjDckDcnl det[J e ] I l„cWlW m Wc^?jk. 
ijk c 

Using the properties of the Kronecker delta while contracting on the indices i, j 
and k gives 

.? AUjkfe= (6-6.12) 

ijk 

£ (r|r| + r£r* + r|rz)amnU£iD£i| det[J e ] I amnWaWmWn^mn 
ia 

+ ? (rfsi + r£s§ + r|s|)ibnI>IiD&«|det[J e ]|ibnWiwbW n ^bn 

ib 

+ ? (r|t| + ryty + rit|)i B cDIiD£ n |det[J e ]|i BC wiw l nWc^? B1 c 
ic 

4 E (Sx^x 4 SyTy 4 sfr®)aiinDBjDal I d6t[J e ] | a»n w a w m w n$ajn 
j * 

4 E (s£s£ 4 4 sIsDlbnD&jDgBldet^llbnWlWbWn^fjn 

j*> 

4 E (Sxt® 4 Syty 4 s|t|)lncI^Bj^cn | dct[J®] | lmcWl^n^<$5jc 
j c 

4 E (tx^x 4 tyrS 4 tfr®)aain^nkl^al| d®^[^ ] | anmWaWnWn$amk 
ka 

+ £ (t£s£ + tSs$ + t^)ibnD£kDEa|det[J e ]|lbnWlW b W n ^bk 
kb 

4 E (txtx 4 tyty 4 t®t®)lmcl^ck^cn | dct[J € ] | ImcWlWniWc^emk* 
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The superscripts on D indicate the local direction associated with the 
derivative. We must maintain this designation, since, in general, the number of 
nodes in each of the three local directions differ. In other words, the quantities 
D r , D s , and D l differ. A simplification occurs upon introduction of the 
symmetric array of coefficients, 

(Cll)aan = (6-6.13) 

(r|r| 4~ ^yly + r|rz)amn|d6t[J e ]|aan^a^«^n> 

(CnWta = (6-6-14) 

(rjs| + + r|s!)ibn I det[J'] | ibnWiwbW„, 

(C 13 )f.c = ( 6 - 615 ) 

(ixtx + fyty + ^ztz)lmc | dct[J e ] | lac^l^n^C) 

(CzOU H ( 6 - 616 ) 

(s|r| + Syly + Szr|)amn | d6t[J € ] | amnWaWmWiii 

(C 22 )?bn = ( 6 ‘ 617 ) 

(8|sl + SySy + s|s|)lbn|det[J e ]|lbnWlW b Wn, 

(C 23 )W = ( 6 - 618 ) 

(sit! + S$t£ + sltDlmc | det[J e ] I lmcWlWaWc, 

(C,0S.» 3 (6-6-19) 

(t|r| + tfrf + t|t|) amn | det[J e ] | amWjWiW n> 

(C 32 )?bn = ( 66 - 2 °) 

(tx^x 4 " tySy 4 “ t|Sz)lbn | det[J e ] | IbnWlWbWm 
and (C 33 )5mc = (6.6.21) 

(txtx ■+■ tyty + tztz)lnc | det [J e ] | IbcWIWbWc, 
into (6.6.12), giving 

.? Afmnijk $ijk = (6.6.22) 

ljk 
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? (Cii)lmn DaiD£i$imn + ? (C| 2 )fbnDIiD§ B |)ibn + 
ia ib 

E (C i3)WDfiD engine + ? (Ci2)a«nD£j®al^ajn + 

ic J a 

? u (C 2 2)?bnD| j D|^! jn + ? (C 2S )f.cDSjDyfoc + 

J 0 JC 

S (Ci3)!nmD nkDal^aak + E (C2 3 )fbnDnkIHm^bk + 

ka kb 

E (C 33 )?-cDkD^f-k. 

kc 

Reordering the terms in this relation yields 

S Af„iikfe= (6.6.23) 

ijk 

S (CnJImnDiiD^U + E (C l2)a«nDalDmj^ljn + 

ia ja 

E (Ci3)|nn DalDnk^amk + E (Ci 2 )fbnD§*DTi$ibn + 
ka ib 

E (C22)!bnD&jD| m $5 jn + E (C 23 )?bnD§ B D‘ k $?bk + 

j b kb 

E (C 13 )fmcD£nDIi$? mc + E (C 2 3 )?mcD^DS j fcc + E (CasJf.cD^D^k- 

ic jc kC 

Further arrangements produce 

.? A? Bnijk $? jk = (6-6.24) 

ijk 

ECDDH? (C n)amn(l-^ai^inin) + ? (C 12 )WDSjfen) + £ (C 13 )i«n(D£ k ^i 1Ilk )] + 
E[E (C 12 )!bn(DIi$? b n) + ? (C 2 2 )?bn(D| j fc„) + E (C 23 )fbn(D£ k $?bk)]D§ B + 

b i J * 

E[E (CuJWDIifcc) + ? (C 23 )? BC (DSjfcc) + ? (C33)?»c(D^«k)]D^.t 

c i J * 

Introducing 

BaSn = ?Bai$iiin> (6.6.25) 

i 

BISn = pjjfen, (6 6.26) 

and BIS. 5 IDjkfc.k, (6.6.27) 

k 

and substituting into (6.6.24) gives 


JThe quantities within parentheses, for example p a i4ian> can be written as A a m 
after performing the matrix multiplication. In this form the product D a iA aBn 
does not satisfy the rules of matrix multiplication. In matrix multiplication the 
second index on D a i must equal the first index of Aam- This involves 
reorienting the indices of D a i by using the transpose. Hence the product should 

be written as (Dia^A^n, where the superscript t refers to the transpose. 
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(6.6.28) 


E Afanijk $tjk — 

ljk 

E(DI a ) t [(Cn)| B .nB£Sn + (C,2)l» n B|£n + (C, 3 )i«nBiSn]+ 

a 

E[(Ci 2 )?bnBIg n + (C 22 )?bnBf!n + (C 23 )!bnBt| n ]D| B + 
b 

E[(C, 3 )! bc BT£c + (C 23 )f„ c B?Sc + (C 3J )f BC Bi£ c ]Dcn- 

C 

This equation applys throughout a single element. Equation (6.6.7) summed 
over all elements gives 

E .? [Afnnijk — A 2 n?nnijk]$ijk = 
e = l ljk 

E E ; n? Bnijk F1 jk . (6.6.29) 

e = l ijk 

This matrix expression for the unknown expansion coefficients $f mn appears in 
the standard form, Ax=B. These coefficients represent the velocity at the 
nodes. We use this expression in a global node numbering system, where the 
global nodes lying on faces shared between adjacent elements receive a single 
label and contain a single value for each variable, though several elements may 
contribute. The assembly process, referred to as ‘direct stiffness’, properly 
accounts for the contributions from several elements. 
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APPENDIX A 

CO-ORDINATE TRANSFORMATION 
A.l Global and Local Co-ordinates 

We decompose the given domain into arbitrarily shaped six-sided elements, 
we then transform these elements from the global co-ordinate system which 
applies to the entire domain to a local co-ordinate system within each element. 
The elements in physical space transform into cubes in local space, with local co- 
ordinates ranging between ± 1. Since the physical equations apply in global 
space, we must develop a method for representing spatial derivatives and 
integrals in the local system. Specifically, we require a method to transfer from a 
co-ordinate system (xi,X 2 ,xs) to a system (£ 1 ,^ 2 , &)• 

Assuming an expression such as 

Xi = Xi(£i,(2>£s) (A.1.1) 

exists between the two co-ordinate systems, the total derivative of Xi appears as 



Expressing this relation as a vector by letting dx* equal the vector with 
components (5xi/$£j)d£j, changes the total differentials into 



dx 1 = 

(A.1.3) 


dx 2 = |^d6 e i + ^d^a + ^d£3e 3 , 

(A.1.4) 

and 

d* 3 = ^d£i e i + 

(A. 1.5) 


We require differentials in vector form since the various surface integrals 
encountered in the analysis contain differential areas pointed in the surface 
normal direction. 

A.2 Transform Between Global and Local Co-ordinates 

Equation (A. 1.2) represents a valid procedure for transforming from the 
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global co-ordinates (x,y,z) to the local co-ordinates (r,s,t). In matrix form this 


transformation equals 


’dx 


X r x s x t ' 


dr' 

-w 

dr 

dy 

= 

y r y 8 yt 


ds 

ds 

dz 


_Z r Zg z 


dt 


dt 


where the Jacobian, [J], equals 



X r X S X t 

y r ysyt • 

Zg 


(A.2.1) 


(A.2.2) 


Subscripts refer to partial derivatives with respect to the local co-ordinates. The 
inverse transform, whereby the local co-ordinates appear in terms of global ones, 


equals 


'dr 


dx' 

ds 

= [j- 1 ] 

dy 

dt 

dz 


The inverse of the Jacobian equals 


l J_1 l = dilTJJ 


yszt-ytzs 

ytzr-yrzt 

.yrZs-ysZr 


which simplifies, upon introduction of 


x t zs-x s zt 

X r Zt-XtZr 
XgZf x r z s 


x s yt-xty s 
xtyr-x r yt , 
x r ys-x s yr 


(A.2.3) 
(A. 2.4) 


(A.2.5) 



xi = ysZt-ytZs, 

x 2 = ytzr-y r z t , 

X3 = yrZs-ysZr, 


yi = x t z s -x s z t , 

y 2 = XrZt-xt z r , 

y3 = XsZf— X r Z s , 

and 

zi = xgyt-xtys, 

z 2 = x t yr-x r yt, 

z 3 = Xrys-Xsyr, 


to 



xi yi zi 

X 2 y 2 Z 2 • 
X3 y3 Z3. 


(A. 2.6) 


Substituting into the inverse transform results in 


dr" 

ds 

dt 


1 

"detpi 


xi yizi 

x 2 y 2 z 2 
Xs y3 Z 3 . 


dx 
dy » 
dz 


(A.2.7) 


which represents a convenient method for computing total derivatives of local co- 


ordinates as f un ctions of global differentials. Partial derivatives of the local co- 


ordinates occur when two of the global co-ordinates remain constant, or, 
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equivalently, by setting one or more of the differential quantities dx, dy, or dz 
equal to zero. 

A.3 Volume Transformation 

The variational form of the equations of motion contain volume integrals in 
physical space; however, the expansion functions only apply in transform space. 
Consequently, we require a transformation procedure for conversions between the 

two systems. 

If the vector dx 1 represents a differential length in the global system, a 
differential volume equals 

dV = dx l »(dx 2 *dx 3 ). (A. 3.1) 

Using the previously derived expressions for the vectors dx 1 gives 

dV = Cijkfjrdfi fjrj^i (A.3.2) 

where tijk equals the alternating unit tensor. Introducing the Jacobian gives 

dV = det[J]d£,d£2d6; ( A - 3 - 3 ) 

where the determinant of the Jacobian depends on 
det[J]d£id6df3 = fijk||j||y 

(A3 ' 4) 

We require 

0 < | det[J] | < m (A. 3. 5) 

for a proper transformation. With (A.3.3) we have a method for transforming 
volume integrals between the global and local co-ordinate systems. This 
transformation appears in the variational statement of the governing differential 
equations, or whenever a volume integral in the global system appears. 

A.4 Partial Derivatives 

The governing differential equations require partial derivatives of the 
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dependent variables in terms of the global co-ordinates. Since the global co- 
ordinates depend on the local co-ordinates, and vice versa, partial derivatives 
may appear as 



dp _ dpdx , dpd 8 , dpd t 
cTx didx ifsitx dtdx’ 

(A. 4.1) 


dp dpd i , dpd s , dipdi 

(A.4.2) 

and 

dy> _ dpdx , dpih , dud t 
dz dxdz + dsdz ^ dtdz' 

(A.4.3) 

The partial derivatives of local co-ordinates depend on (A.2.7). 
the partial derivative dr/dx equals 

For example, 


d r dr 1 _ 

EE = EE y,z-3etpj Xl ' 

(A.4.4) 

Using similar expressions for the other derivatives gives 



dip _ xi dtp , X 2 dip , xj dip 
det[JJ0r ^ detlJJfls ' detlJjdt 1 

(A.4.5) 


d<P- 7i du , y 2 dip j % du 
<Jy det[Jj5r ' det [ Jjds 1 det[Jjdt’ 

(A.4.6) 

and 

dip z\ dp , z 2 dp , z% dp 

M ~ det[JJ?r + det[JJ3s + det[j]3f 

(A. 4.7) 


Expressions such as these appear whenever we represent partial derivatives of the 
dependent variables in transform space. 

A.5 Partial Derivatives Using Spectral Expansions 


The approximation procedure expresses the dependent variables as a 
truncated series of interpolating functions and coefficients. Therefore, 
expressions such as (A. 4. 5), containing the dependent variable p, actually 
contain a series of coefficients and functions which depend on the local co- 
ordinates. Our method uses interpolating functions based on the Legendre 
polynomials. For example, the one-dimensional interpolation function of order n 
for the dependent variable p equals 

I a ip =l<pkM T )- (A.5.1) 

k = 0 

which contains the interpolating functions, {^(r), k=0,l,...,n}, and coefficients, 
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fa. The derivative of <p with respect to the local variable r appears as 

DnvK 1 ) = k ^ 0 Wcgfl^( r )]- (A.5.2) 

In the current formulation, we compute the derivatives at the nodes and not 
between nodes. Evaluating the derivative at r a gives 

D n y<ra) = k S o v=>kDak, (A.5.3) 

where gj{^k(r)] r =a=Dak represents the derivative of the expansion function ^k 
with respect to the local co-ordinate r at node r a . Substituting similar 
expressions for the other partial derivatives in (A. 4.5) gives 

(|f)Ibc = (det[JJ^ bc i ? k ^i k ® ai ^!i( Sb )^ k ( tc ) (A.5.4) 

+ ( ^et[Jj ) abc •? k ^ijkV'i( r a)Dbj V’k(tc) + (dei[Jj) abc i ? k ¥ , ij k V’i( r a)^j( s b)Dck- 
The superscript e represents the element number, while the subscripts abc 
stand for the spatial node under consideration. Since the interpolating 
polynomials satisfy ^i(r a )=^ai) the partial derivative expressions reduce to 

(G.5.5) 

(|f)Ibc= ( det[JJ ^ bc ? ^ bcDai + Uet{jp abc f ^J' cDb i + (dei['JI^ bc k ^ bkDck ’ 
(|^)lbc= ( de t[Jj)* bc i ^ bcDai + (detfjj^^ f ^ jcDbj + Qet[j|^ bc k ^ abkDck ’ 

(|f)Ib C = (Heffjj)^ bc ? ^ bcDai + (def[JJ^ bc f ^ jcDbj + (de?pj)* bc k ^ bkDck ' 

These equations represent the partial derivatives of the dependent variable, <p, 
at the nodal points. At this point we impose no restriction on the location of a 

node within an element, i.e., the node may lie in the interior or on the surface of 

an element. The next section computes derivatives along element faces by 
applying (A.5.5) to nodes lying along the surface. 

A.6 Derivatives Along Element Faces 

The calculation of surface stresses requires derivatives along element faces. 
The equations for derivatives at any node (a,b,c) simplify along the element 
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surface since one of the indices remains fixed while the other two vary. This 
results in a set of expressions for each of the six dement faces: 

(^)tbi= (aerjjj)tbi s £?biDai + (gg£pj)*bi ? ^ji D bj 

+ (a^fjj)Ibi £ vlbkDjk, (A.6.1) 

(|^)Ibi= (<Jetjjp abl ? ^ blDai + (aet[Jl)* bl ? ^ii Db i 
+ (det[Jp abl k ^ bkDlk 

(|f)abl= (g^(jj)abl S ViblDai + (g^pj)Sbl ? ¥>tjl D bj 
+ (^jj)IbiS ^IbkDik, along the side at t= — 1; 

(|^)Ijtc= (det[jj)^' tc ? ^J' tcDai + Qlet[Jp aitc f ^ cD j tJ ' 

+ (det [Jp ajtc k ^J tk ^ ck (A. 6. 2) 

(|y)ajtc= (det[Jj) ajtc ? ^ tcDai + (det[Jp ajtc f ^ ajcDjt j 
+ (det'[Jp a -> tc k ^J‘ tkDck 

(|^)ajtc= (difpj)ajtc £ ^fjtcDai + (^fjjJajtc ? VajcDjtj 

+ ( de ?lJj)ai tc jj* ^ajtkDck) along the side at s=l; 

(|^)abkt= (diT[J[)abkt £ ^bktDai + (g^pj)abkt ? PajktDbj 

+ (def['J]) abkt k ^ bk ® ktk (A. 6. 3) 

(|^)abkt= (ditjjj)abkt £ ^bktDai + (det[Jp abkt f ^i ktDb J 
+ ^det [Jp abkt k ^ abkDktk 

(§£)abkt= (deffjj) abkt i ^ bktDai + (a^[3I )Sbkt f V’ajlaDbj 
+ (g^jjj)abkt ^ JbabkDktkj 3long the side t=l; 

(^)aic= (det[Jp alc ? ^ lcDai + (def[Jp alc J ^ ajcDlj 

+ (det[jj )Ilc k ^ alk ^ ck (A. 6.4) 

(|^)aic= (deffjj) alc ? ^ lcDai + ^det[Jp alc f ^ a i cD >i 
+ Qlet[Jp alc ? ^ alkDck 

(|f)aic= (deljjpaic ? V’ticDai + (ge?fjp alc f ^ajcDij 
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(A.6.5) 


+ HeffjJ^ lc ? ^ lkDck ’ the siie 8 = “ v » 

(|f)!l>c= (a^[ 3 j)?bc s V’tbcDii + (g^fjj)!bc ? ^jcDbj 

+ (3eT[JJ^ bc k ^ bkDck 

(|^)1bc= (aifjjj)!bc S ^bcDii + ( 3 ^jjj)?bc S £fj c Dbj 

+ (3etfjJ^ bc k ^ bkDck 

(|f)1bc= (a^jjj)tbc S V’tbcDii + (gefpj)fbc j ^j C D b j 
4* ( ^ e ^j] )1bc S ^Ibk^ck) along the side r= — 1, and 
(|^)?tbc= (aefjjj)* tbc ? ^bcDiti + (^[Jj)itbc s ^itjcDbj 

+ (a^fjjKtbc s ^tbkDck (A. 6. 6) 

(|^)itbc= (3etfjj) ?tbc ? ^bcDiti + ( 3 ^pj)itbc £ ^tjcDbj 

+ (det[j]^ tbc k ^ tbkDck 

(|f)ttbc= (aHpj)*tbc s ^bcDiti + (aifpj)*tbc ? ^tjcDbj 
+ (^jj)itbc s £!tbkDck, along the side r=l. 

A.7 Surface Integrals 

Often surface integrals involving the differential area vector appear in 
boundary conditions and computations of surface stresses inter alia. According 
to (A. 1.3) through (A. 1.5), the transformation between global co-ordinates 
(x,y,z) and local co-ordinates (r,s,t) equals 

dx = ^dre r + ^^ Se s + (A.7.1) 

dy = |^dr er + Ifdses "** it^ tet * (A. 7.2) 

and dz = ^dre r + ^dses + ^dtet- (A. 7. 3) 

In the local co-ordinate system, the element surface depends on two of the three 
local variables while the third remains fixed. For example, along side one co- 
ordinates r and s vary, while co-ordinate t equals —1. Hence, the 
transformation along side one appears as 
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dx 1 = ^jd re r + ^g^ se ^> (A. 7.4) 

dy 1 = + fjrdses. (A.7.5) 

and dn^^drer + ^dseh, (A.7.6) 

where here the superscript refers to side one and not direction one as in (A. 1.3). 
After integrating these relations, expressions for the physical co-ordinates in 
terms of the local co-ordinates along the surface appear: 

xi = x(r,s), (A.7.7) 

y 1 = y(r,s), (A. 7.8) 

and z 1 = z(r,s). (A.7.9) 


Let R^x^x+yiey+ziez represent a position vector describing the location of 
any point along surface one. Substituting the previous expressions for the 
physical co-ordinates into this position vector gives 

R 1 = x(r,s)e x + y(r,s)e y + z(r,s)ez. (A.7.10) 

A differential area vector in the normal direction to the surface corresponding to 
a differential change in the local co-ordinates r and s equals 

dA 1 = — (A. 7. 11) 


The minus sign indicates that the area vector points in the local coordinate 
direction — et- We can write the two partial derivatives as 



dR dRdx„ , dRdy„ . dRdz„ 
W = dxdi** '* + bzfc**' 

(A.7.12) 

and 

dR 5Rftc„ , dRdy„ , dRdz^ 
JT = W8 s** + + 

(A. 7.13) 

Using (A.7.10) reduces these expressions to 



5R 

jY — Xre x + yie y + ZjCz, 

(A.7.14) 

and 

3R 

jY = Xse x + yse y + z^. 

(A.7.15) 


Substituting these expressions into dA 1 and taking the cross-product gives 
dA 1 = - [(y r z s - z r y s )e x + (z r x s - x r z s )e y + 
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(x r y s - yrx s )ez]drds, (A.7.16) 

or, using the simplified notation of (A.2.5) 

dA 1 = — (x 3 e x 4- J&y + * 3 e®)drds. (A.7.17) 

Similar expressions for sides two through six equal 

dA 3 = (x2e x + y2ey + z^tlrdt, (A.7.18) 

dA 3 = (x 3 e x + y 3 ey + z^drds, (A. 7. 19) 

dA 4 = — (x2e x + y2e y + z^Jdrdt, (A.7.20) 

dA 5 = — (xie x + yiCy + zie^Jdsdt, (A. 7. 21) 

and dA® = (xje x + yi6y + zje^Jdsdt. (A. 7.22) 


Expressions such as these will appear whenever we compute surface integrals. 
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appendix b 

SURFACE FORCE AND MOMENT CALCULATION 
B.l Background 

Fluid flowing past solid boundaries develop, surface stresses winch lead to 
surface forces and moments after integration over the surface area. The elements 
located adjacent to these surfaces contain flow-field information necessary for 

surface stress calculations let dA equrf a differential area on the element 

surface with unit normal vector m Let the force per unit area due to the flmd 

stresses on this surface equal t(n). Then, the total force on a finite surface, M. 

results from integrating the contributions from all the differential elements: 

F=/f t(n)dA. (BU) 

Similarly, the total moment about a given point due to the surface forces equals 

M = // X*t(u)dA, (B12) 

JJ dn 

where x represents the distance between the surface element and the point 
about which the moment is taken. Let 0 ji denote the i th component of t(j), 

and let t(n)i denote the i th component of t(n). Whence 

, % (B.l. 3) 

t(n)i = 0 ji n j> v 

where the stress tensor equals 

an a 12^13 (B.1.4) 

ff = (721 022 ^23 • 

~ 031 032 

Since t(n) represents a vector and n represents a unit vector independent of 

the 0 ji, we require 

t(n) = n-£ (B.1.5) 

This equation represents the force per unit area due to fluid stresses on a surface 
with unit normal vector n. When this expression replaces t(n) in the force and 
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moment integrals, terms such as n-odA and arise. Since dA 

represents a scalar, we may write these expressions as ndA* a and x*(ndA*<r) 

We developed expressions for ndA=dA along each of the six element faces in 
appendix A. We still need, however, an expression for the stress tensor. 

B.2 Incompressible Stress Tensor 

In Cartesian co-ordinates the diagonal components of the incompressible 


flow stress tensor equal 



a xx = — p + a yy — — P + ^cly’ 


and 

far 

Ozz = — P + 2/^-, 

(B.2.1) 

while the off diagonal terms equal 



,d u , dv\ 

a X y = <7yx = IVqfj + 

(B.2.2) 


id u , dw>. 

Vxz = 0 Z x = A-fo + -fch 

(B.2. 3) 

and 

f dv , dvrs 

<Jyz — °zy — + fjj)- 

(B.2.4) 

B.3 Surface Force 



As indicated in 

§B.l, the force over a surface dCl due to 

the fluid stresses 

equals 

F = ff t(n)dA= ff n-udA. 
JJ d(l JJ dCl ~ 

(B.3.1) 

The components of the force in the global (x,y,z) system equal 



F x = ff (<r xx n x + (Tyx^y + <7zxn z )dA, 

JJ d n 

(B.3.2) 


Fy = ff ( CT X yn x + CTyyUy + <7zyUz)dA, 

JJ dn 

(B.3.3) 

and 

F z = ff (<r X z n x + ^yzUy + <r zz n z )dA. 

JJ d n 

(B.3.4) 

Since n x dA, nydA, 

and n z dA represent the projection of dA 

onto the planes 


perpendicular to the global co-ordinate axes, we may write n x dA— e x *dA, 
nydA=e y -dA and n z dA=e z -dA. Substituting these expressions into the force 
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components gives 



F x — ff (<^xx^x + GyxGy + 0zxGz)"dA, 

(B.3.5) 


Fy = JJ* (<Txy^x + <7yyGy + ^yCz)*dA, 

(B.3.6) 

and 

F z = ff (cTxzCx + ^yzCy + 0zz£z) *dA. 

•'•'an 

(B.3.7) 

The forces along side one, where dA l = — (x 3 e I +y 3 e y +Z 3 e z )drds, equal 



Fx 1 — — /*/* ( Oxx x 3 4" 0yxy 3 + ^zx z 3 )drds. 

JJ dw 

(B.3.8) 


Fy 1 = - ff (a xy X3 + ffyyjz + <r Z yZ 3 )drds, 

(B.3.9) 

and 

Fz 1 = — ff (a xz x 3 + a yz y 3 ■+■ 0 zzZ 3 )drds. 

JJ dW 

(B.3.10) 


On the local level, the surface of integration reduces to a square with sides 
ranging between ±1. Using the spectral collocation approximation allows an 


expression for the force Fx 1 as 

Fx ie = f \f V**x 3 + + 

1J " -1 

o r zxZ3)iji^i(r)V'j(s)drds > (B.3.11) 

where the expression (...)iji indicates that the terms within the parentheses 
depend on the local node yi of element e. Integrating over the element surface 
gives 

F x ie = — E (OxxX3 + Oyxji + tfzxZsJijiWiWj- (B.3.12) 

ij 

In a similar manner, the y and z components appear as 

F y l e = -E (<r xy X 3 + (T yyY 3 + tf 2y Z 3 )?j iWiWj (B.3.13) 

and F z ie = — E. (oxzX 3 + 0yzY3 "b 0 zzZ 3 )ijiWiWj. (B.3.14) 

The forces on faces two through six equal 

Fx 2e = E^ (o"xx x 2 + ^yxy 2 + <7zx z 2) ljtkWiWk, 

Fy 2e — E^ (CxyXj + CTyyY 2 + <7zy Z 2) lj tk^iWlt, 

and F z 2e = E^ (ctxzX 2 + 0 y Z y 2 + 0 ’zzZ 2 )ijtkWiWj ( , (B.3.15) 
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on face two, where dA^x^x+y^y+z^^dt; 

Fx 36 = ?. (^xxX3 + 0yxy3 + Gzx z 3)ijkt w i w j> 

Fy 3€ = E (<TxyX3 + GyySl + 0 fcy z 3 )ijktWiWj, 

ij 

and F z 3 e = E (a xz x 3 + + ^sJtjktWiWj, (B.3.16) 

ij 

on face three, where dA=(x 3 e x +y 3 fcy+ z 3 Cz)drds; 

F x * e = E (^xx^2 + ^yxj2 + 0zx z 2)ilkWiWk, 

ik 

Fy 4e = E (axyX 2 + ayyy 2 + a zy z 2 )?ikWiWk, 
j ik 

and F z < e = E (<r xz x 2 + ^2 + ^2)1ikWiW k> (B.3.17) 

ik 

on face four, where dA= — (x 2 e x +y 2 e y +Z 2 e z )drdt; 

F x 5e = E (a xx x t + ayotfi + a zx zi)!jkWjWk, 
j* 

Fy 5e = E (cr xy Xi + ffyyj 1 + Oky z l)tjk w j w k> 
jk 

and F z s e = E (a xz X! + ayzyi + azzzOljkWjWk, (B.3.18) 

jk 

on face five, where dA= — (xjei+yiey+zie^dsdt; and 

F x fle = E O xx X! + (Tyx yi + azxZl)itjkWjWk, 
jk 

Fy® e = E ( O x yX 1 + (7y y y 1 + <7zyZ l) 1 tj k w j W k , 
jk 

and F z 8 e = E (a xz xi+ cr y zyi+ ^zzZi)?tjkWjWk, (B.3.19) 

jk 

on face six, where dA=(xie x +yie y +zie z )dsdt. Note that when any of the 
subscripts ijk equals 1, we refer to the node at the face corresponding to rst 
equal to —1, and when any of the ijk equals it, jt, or kt, we refer to the 
node at the face corresponding to rst equal to 1. Also, the basis functions 
integrated over the local co-ordinate domain, wi= j*^i(r)dr, wj= /^j(s)ds, 
and Wk = f Vi(t)dt, differ when the number of nodes in the r, s, or t 
directions differ; consequently, we must remember that Wi refers to the r- 
direction, Wj refers to the s-direction, and w k refers to the t-direction. 

B.4 Surface Moment 
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The moment along a surface due to the fluid stresses appears as 


M = ff y*n*adA, 
JJ dQ 


(B.4.1) 

where the moment "arm" equals ^XxCx+Xyfy+X*®!- Substituting into (B.4.1) 


and performing the cross product yields the moment components in the global co- 
ordinates as (B 4 2) 

Mx — ff [Xy( ^XzUx + Pyzlly + (TzzPz) “ Xz(^xy n x + O'yyRy + <7 Z yn z )]dA, 

"dtl 

= SL [Xz( cr xxHx + 0yx% + 0zxBz) “ XxC^xz^x + ^yz&y + ^zz^ 2 )]dA, 
dQ 

M z = ff [XxC^xy^x + ^yy n y 4 0zy n z) Xy(^xx n x “I" GyxHy 4 Gzxn z )]dA. 
"dQ 

Using the same notation as in the surface force calculation gives 

M ‘ = XL [Xy(^xzCx 4 (TyzBy 4 0zz®z) Xz( ^xyCx 4 0yy€y 4 0zyCz)]*dA, 
dQ 

My = ff [Xz(^xxCx 4 ^yx€y 4 (Tzjfi z) “ Xx(^xzCx 4 ^*yz®y + 0zzCz)]*dA, 
"dQ 

M ‘=/C, [Xx(^xy^x "I" 0yy®y 4 0zyCz) Xy(^xx£x 4" Gyx^y 4 (7 Z x€z)]"dA. 

* dQ 

(B.4.3) 


For example, the x-component evaluated along side one, where 
dA=— [x 3 ex+y 3 ey+z 3 ez]drds, equals 

Mx 1 ® = — J j[Xy( GxzX-i + ^yzy 3 + o'zz^s) ~ Xz(&xy x 3 + &yyY 3 + cr zy z 3 )] e drds. 

Introducing the spectral expansion functions gives 

Mx ie = —E J" ^J“ |[Xy( <7 xzX3 + (Jyzy 3 (B.4.4) 

+ OzzZ 3 ) — Xz(<^xy x 3 + <7yyy 3 + ^zy z 3)] lj ( s) dxcis . 

After integration we obtain 

M X 1C = -S [Xy(ffxzX 3 + Oyzy 3 + ffzzZ 3 ) 

- Xz(^xyX 3 + a yy y 3 + tr zy z 3 )]ijiWiWj. (B.4.5) 

In a similar manner, the y and z components appear as 
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My ie = — E [Xz(^xx x 3 + O’yxy 3 + &zx z z) ~ Xx(<7 X zXs + VytJS + <7zzZ3)]!jlWiWj, 
ij 

M z * e = — E [Xx(tfxy x 3 + &yyY Z + Ozy z 3) — Xy( ff xx x 3 + <*yxJz + ^zx z 3)]ij l w i w j- 
ij 

The moments along faces two through six equal 

M x ae = E^ [Xy(<7xzX2 + 0yzy2 + £7zzZ2) “ Xz(^xy x 2 + VyyYl + ^zy z 2)]ijtkWiWk, 

My je = E [Xz(^xx x 2 + &yxY 2 + & zx z 2 ) — Xx(^xz*2 + OyiZl + ^zz^jijtkWiWk, 

ik 

= E [Xx(o’xy x 2 + Oyyj2 + ^zy z 2) “ Xy( ff xx*2 + Oyrjl + 0 r zx z 2 )]ijtkWiWk, 
ik 

along face two, where dA=(x 2 e x +y 2 e y +Z 2 e z )drdt; (B.4.6) 

M x 3e = E. [Xy(^xzX3 + QyzY Z + 0zz z 3) ~ Xz(^xy x 3 + OyyYz + <7 zy z 3 )]fjktWiWj, 

M v 3e = E [Xz(^xx x 3 + VyxYZ + OzxZ3) — Xx( a xzX3 + Oyi3 Z + OzzZ3)]ijktWiWj, 

ij 

M z 3e = E [Xx(^xyX3 + GyyYZ + 0"zy z 3) “ Xy( ff xx x 3 + a yxYZ + ffzx z 3 )]fjktWiWj, 

along face three, where dA=(x 3 e x +y 3 e y +Z 3 e z )drds; (B.4.7) 

M x 4e = — E [Xy( cr xz x 2 + OyzY 2 + 0 r zzZ2) — Xz(^xy*2 + &yyY2 + 0 'zy z 2 )]ilkWiWk, 
ik 

My< e = -E [Xz(<^xxX2 + OyxY2 + O lx Z<i) - Xx(^xzX2 + OyiSl + ^zz z 2 )]?lkWiWk, 

J ik 

M z 4e = — E [Xx(flxyX2 + 0yyy2 + 0zy z 2) — Xy(°xx x 2 + &yxj2 + ffzx z 2 )]?lkWiWk, 
ik 

along face four, where dA= — (x 2 e x +y 2 ®y+ z 2 ©z)drdt; (B.4.8) 

M x 5 e = -? k [XyC^xzXi + ffyzyi + <7 zz Z 1 ) - Xz(^xyXi + ^ yy yi + <T Z y z l)]1jkWjW k , 

My5 e = -E [xz(^xxXi + a yx yi + a zx Zi) - Xx(^xzXi + ^yzyi + ^zz z i)]tjkWjW k , 
jk 

M z 5e = — E [Xx(tfxyXi + CTyyYl + <7zy z l) _ Xy^xx^l + VyxY 1 + 0zxZi)]fjkWjWk, 

jit 

along face five, where dA= — (xie x +yie y +zie 2 )dsdt; (B.4.9) 

M x fle = E [Xy(^xzXi + (Tyz Yl + <7zzZl) - Xz(^xyXi + ffyyyi + O'zyZi)] itjkWjWk, 

jk 

M y fle = E [xz(<7xxXi + ffjTtyi + a zx zi) — Xx(<^xzXi + ffyzyi + ^zz z i)]itjkWjWk, 

jk 

M z fle = E [Xx(flxyXi + (7 yy yi + CT zy Zi) — XyC^xxXi + CFyxy 1 + 0zxZl)]!tjkWjWk, 

jk 

along face six, where dA=(xie x +yie y +zie z )dsdt. (B.4. 10) 
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APPENDIX C 

IMPERMEABLE WALL BOUNDARY CONDITIONS 
C.l Physical Basis 

Along an impermeable wall the normal fluid velocity must equal the normal 
component of the wall velocity, while the tangential velocity may vary. We 
express this as 

n-V = n-V wa ii. (€.1.1) 

Another condition on impermeable wall boundaries requires zero tangential 
stress, or, equivalently, that the stress points normal to the surface. In other 
words, the cross product of the surface unit normal and the local stress vector 
equals zero. Thus, an impermeable wall, unlike a solid wall, cannot support a 
tangential stress. In symbolic form this boundary condition equals 

n*(£*n) = 0, (C.l. 2) 

where n represents the surface outward unit normal vector, and a represents 
the total stress tensor consisting of the sum of the isotropic and deviatoric terms. 
We must prove that the combination of these two boundary conditions equals the 
single condition given by 

(n-V)V = 0. (C.2.1) 

The derivation of this relation follows. 

C.2 Mathematical Derivation 

Begin by expressing the unit normal as 

n = n x e x H - Uyey H- n x Cz (C.2. 2) 

and the velocity as 

V = Ue x + Ve y + We*. (C.2.3) 

The inner product n-V gives 
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n-V = Un x 4- Vn y + Wn z . 


(C.2.4) 


Setting this equal to the normal velocity of the boundary gives 

Un x + Vn y 4- Wn z = n-Vwaii- (C.2.5) 

The stress vector on a surface with normal, n, equals t(n) which appears as 

t(n) = £*n = (C.2.6) 

^yx% -H 0z X n z )e x + 

( cr xy n x + ayyn y *4“ a zy n z )e y 4" 

(c xz n x 4" ffyjDy 4" ctzzUz)Gz 

according to appendix B. The cross product of this relation with the unit normal 
equals 

n*t(n) = nx<r*n = (C.2.7) 

[n y (c xz n x 4" <r yz n y 4- (Jzz&z) — n z (cr xy n x 4- °yy n y ■+■ <r zy n z )]e x 4- 
[u z (< 7 xx n x 4- 0 y x n y 4- cr zx n z ) — n x (cr xz n x 4- a yz n y 4- cr zz n z )]e y 4- 
[n x ((T xy n x 4- er yy n y 4- a zy n z ) — n y (cr xx n x 4- <Xyxn y 4- c zx n z )]e z . 

Using (B.2.1) through (B.2.4) for the stress components gives the following 
relations for the three components of the above expression: 

(n*£-n) x = /in x n y (|j4-|^) 4- ^y n y(jy+j^) + ^ n y( 2 ^Y~P) ~ 

/xn x n z (|y4-|Y) “ ^y^Jf ~ P) ~ (C.2.8) 

(n*£-n) y = /m x n z (2^— p) 4- / m y n z(|y+^) + Z®* 11 z(/ x 1 0 Z ) — 

/m x n x ( ^ 4 - | ^ ) - zm y n x (|y4.|j) " ” p )> (C-2.9) 

(n*£-n) z = zm x n x (^y + ^) 4- Z m y n x(2^y- — p) + Z m ® n x(^y _ *"7z’) — 

~ p) ” / m y n y(f^ + 3y) ~ ^ n y(fl' l 'fy)- (C.2.10) 
The pressure terms cancel in each of these expressions giving 

(nx£-n) x = /m x n y (|^ 4- jj-) 4- /m y n y (|y + j^) + ^zn y (2jj- - 2^j) - 

/mxn z (|y + fc) ~ ItoiPzijz + Jj)> (C.2.11) 
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(nxa-n) y = /m x n z (2|^ - 2^) + /m y n.(|j + j^r) + + |j) - 

^xn x (|| + ) - |m y n x (|^ + £), (C.2.12) 

(n*£-n) z = /m x n x (|y + 1^) + / in y n x(2|j - 2^) + /®*n x (|y + |j) - 

/ m y n y(|^ + §y) ~ + fy)- (C.2.13) 

After further rearrangement we obtain 

(nx£-n) x = -^n z (nr^y + nyyy + n^y-) + + n y^ + n^-J 

+ /<(tt«| f + ^ + »4)(% w - , >» v ). (C.2.14) 

(n«ff*n) y = -fm x ( n x|j + D yfj + n ^f^) + + n z|^") 

+M“4f+-4+ i 4)(».o-"« w ), ( c - 2is ) 

, , , «u , sv . _ «w> , «u . „ sv , „ jw> 

(nx£-n) z = -/znytnx^ + ny^- + n z ^-) + im^nx^y + ny^y + n^-y) 

+ K n *|j + n y|y + n *^)( nx ^ — n yU). (C.2.16) 

The partial derivatives of the unit normal vector with respect to the spatial co- 
ordinates differ from zero on a curved surface. Hence, we cannot interchange the 
derivative operator in the first two terms of each expression with the components 
of the normal vector (e.g., n x |y +n y ^ H n z ^^ f ^(n x U + n y V + n z W)). 
However, the derivative may move outside the components of the normal vector 
if we only consider surfaces with zero curvature. When we restrict ourselves to 
this case, the above relations may appear as 

(nx£-n) x = — /in z ^y(n x U + n y V + n z W) + /my^(n x U + n y V + nzW) 

+ ^ n >4 + n 4^ + n 4>( n * w - ,i ' v )' ( c - 2n ) 

(n*£-n) y = — ^mj^n x V + nyV + n z W) + + HyV + n z W) 

+ A n^ + ny^ + n4)(n*U - n x W), (C.2.18) 

(n«£-n) z = — /niy^(n x U + n y V + n z W) + /m x ^(n x U + n y V + n z W) 

+ faJfe ny TJy nz iz)( nx ^ — n y^)- (C.2.19) 

fl Q O 

Introducing n*V=n x U+n y V+n z W and n *^ =n x^+ n y^+ n z^ into these 
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expressions yields 

(n*£.n) x = Kny|— n 2 $y)n-V + ,m-V(n y W-n z V), (C.2.20) 

(n*£*n) y = —fi{ n xg z ~ n z^) n * V + jm*V(n z U— n x W), (C.2.21) 

(n*£-n) z = _ M n y^f n ^) n '^ + /m*V(n x V-n y U). (C.2.22) 

After introducing the cross product, the three components of (n*£-n) reduce to 

a single vector expression given by 

(n*£-n) = /x[(nxV)n»V + (n*V)n*V]. (C.2.23) 

This equation also appears as 

(n*£*n) = ^n*[V(n-V) + (n*V)V]. (C.2.24) 

Introducing n-V=n-V wa ii gives 

(nx£-n) = /m*[V(n* V w aii) + (n*V)V]. (C.2.25) 


For a general curved boundary moving with velocity V W aii, the quantity 
V(n* Vwaii), does not equal zero. However, since we are limiting our 
consideration to flat boundaries, this term does equal zero. In this case the 

boundary condition reduces to 

(n*£*n) = /m*[(n*V)V]. (C.2.26) 

Since the original boundary condition, n*£*n, equals zero, we require 

/m*[(n-V)V] = p(*-V)(nxV) = 0. (C.2.27) 

This equation is satisfied when the unit normal lies parallel to the velocity vector 
giving n«V=0. This seems like a most unusual condition and we assume that it 
does not hold for most boundaries. The other possibility requires 

(n*7)V = 0, (C.2.28) 

and represents the desired form of the boundary conditions given in §C.l. This 
proves that the two impermeable wall boundary conditions reduce to 
homogeneous natural conditions when we enforce (1) zero surface curvature to 
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ensure V(n« V wa ii)=0, and (2) a surface with unit normal not parallel to the 
local velocity. 
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